def test_testUfuncs1(self):
     # Test various functions such as sin, cos.
     (x, y, a10, m1, m2, xm, ym, z, zm, xf, s) = self.d
     assert_(eq(np.cos(x), cos(xm)))
     assert_(eq(np.cosh(x), cosh(xm)))
     assert_(eq(np.sin(x), sin(xm)))
     assert_(eq(np.sinh(x), sinh(xm)))
     assert_(eq(np.tan(x), tan(xm)))
     assert_(eq(np.tanh(x), tanh(xm)))
     with np.errstate(divide='ignore', invalid='ignore'):
         assert_(eq(np.sqrt(abs(x)), sqrt(xm)))
         assert_(eq(np.log(abs(x)), log(xm)))
         assert_(eq(np.log10(abs(x)), log10(xm)))
     assert_(eq(np.exp(x), exp(xm)))
     assert_(eq(np.arcsin(z), arcsin(zm)))
     assert_(eq(np.arccos(z), arccos(zm)))
     assert_(eq(np.arctan(z), arctan(zm)))
     assert_(eq(np.arctan2(x, y), arctan2(xm, ym)))
     assert_(eq(np.absolute(x), absolute(xm)))
     assert_(eq(np.equal(x, y), equal(xm, ym)))
     assert_(eq(np.not_equal(x, y), not_equal(xm, ym)))
     assert_(eq(np.less(x, y), less(xm, ym)))
     assert_(eq(np.greater(x, y), greater(xm, ym)))
     assert_(eq(np.less_equal(x, y), less_equal(xm, ym)))
     assert_(eq(np.greater_equal(x, y), greater_equal(xm, ym)))
     assert_(eq(np.conjugate(x), conjugate(xm)))
     assert_(eq(np.concatenate((x, y)), concatenate((xm, ym))))
     assert_(eq(np.concatenate((x, y)), concatenate((x, y))))
     assert_(eq(np.concatenate((x, y)), concatenate((xm, y))))
     assert_(eq(np.concatenate((x, y, x)), concatenate((x, ym, x))))

    def concatenate(self,value,axis=0):
        """ Concatentate UncertContainer value to self.
            Assumes that if dimensions of self and value do not match, to 
            add a np.newaxis along axis of value
        """

        if isinstance(value,UncertContainer):
            if value.vals.ndim == self.vals.ndim:
                vals = value.vals
                dmin = value.dmin
                dmax = value.dmax
                wt = value.wt
                uncert = value.uncert
                mask = value.mask
            elif (value.vals.ndim + 1) == self.vals.ndim:
                vals =  ma.expand_dims(value.vals,axis)
                dmin =  ma.expand_dims(value.dmin,axis)
                dmax =  ma.expand_dims(value.dmax,axis)
                wt =  ma.expand_dims(value.wt,axis)
                uncert =  ma.expand_dims(value.uncert,axis)
                mask =  np.expand_dims(value.mask,axis)
            else:
                raise ValueError('Could not propery match dimensionality')
                
            self.vals = ma.concatenate((self.vals,vals),axis=axis)
            self.dmin = ma.concatenate((self.dmin,dmin),axis=axis)
            self.dmax = ma.concatenate((self.dmax,dmax),axis=axis)
            self.wt = ma.concatenate((self.wt,wt),axis=axis)
            self.uncert = ma.concatenate((self.uncert,uncert),axis=axis)
            
            self.mask = np.concatenate((self.mask,mask),axis=axis)
        else:
            raise ValueError('Can only concatenate with an UncertContainer object')

 def computeStatistics(self):
     dates = []
     values = []
     for ts in self.itervalues():
         dates.append(ts.dates)
         values.append(ts.values)
     dates = concatenate(dates).astype(int)
     values = ma.concatenate(values)
     minDate = dates.min()
     nDates = dates.max() - minDate + 1 
     hist = histogram(dates*(-values.mask), bins=arange(nDates)+minDate)
     values = values.compressed()
     if values.size == 0:
         raise Exception("empty dataset")
     lowend = values.min()*0.99
     highend = values.max()*1.01
     #cumfreqs, lowlim, binsize, extrapoints = cumfreq(values, 40, (lowend, highend))
     #normcumfreqs = cumfreqs/values.size
     #ind = ((normcumfreqs > 0.02) & (normcumfreqs < 0.98)).nonzero()[0]
     #if ind.size == 0:
     #    raise Exception("empty dataset")
     #min = ind[0]*binsize + lowlim
     #max = ind[-1]*binsize + lowlim
     #return min, max, hist
     return values.min(),values.max(),hist

def rect2sphere(vector, degree=True):
    """\
    Convert vector (x,y,z) from rect to sphere coordinates. If degree is 
    ``True`` the unit will be in degree.

    Examples
    --------
    >>> convert.rect2sphere([1,1,1], degree=False)
    array([ 0.78539816,  0.61547971,  1.73205081])

    >>> convert.rect2sphere(numpy.array([[1,2],[1,0],[1,3]]), degree=True) 
    array([[ 45.        ,   0.        ],
           [ 35.26438968,  56.30993247],
           [  1.73205081,   3.60555128]])

    """
    x, y, z = vector

    r = np.sqrt(x**2 + y**2 + z**2)
    lon = np.arctan2(y,x)
    lat = np.arcsin(z/r)

    if degree:
        lon = np.rad2deg(lon)
        lat = np.rad2deg(lat)

    return ma.concatenate([ lon[np.newaxis], lat[np.newaxis], r[np.newaxis] ])

    def add_solver(self, fname):

        # Reg exp: Any line starting (ignoring white-space)
        # starting with a comment character. Also col sep.
        comment = re.compile(r'[\s]*[%#]')
        column  = re.compile(self.sep)

        # Grab the column from the file
        metrics = []
        file = open(fname, 'r')
        for line in file.readlines():
            if not comment.match(line):
                line = line.strip()
                cols = column.split( line )
                data = atof(cols[opts.datacol - 1])
		metrics.append(data)
        file.close()

        if self.metric is not None:
	    self.metric = ma.concatenate((self.metric, [metrics]))	    
        else: 
            self.metric = ma.array([metrics])

        # Current num of probs grabbed
        nprobs = len(metrics)
        if not self.nprobs: self.nprobs = nprobs
        elif self.nprobs <> nprobs:
            commandline_error("All files must have same num of problems.")

def fix_pop_grid(tlon,tlat,data):
    """
    Pad coordinates and data on CCSM/POP gx3v5 grid so it can be
    plotted with matplotlib/basemap
    tlon,tlat,data must be 2D arrays

    Inputs:
        tlon, tlat, data

    Outputs:
        lon, lat, data

    """
    # make lon monotonic and pad coordinate and data arrays along lon axis
    tlon = N.where(N.greater_equal(tlon,tlon[:,0].min()),tlon-360,tlon)
    lon  = N.concatenate((tlon,tlon+360,tlon+2*360),1)
    lat  = N.concatenate((tlat,tlat,tlat),1)
    if hasattr(data,'mask'):
        data = MA.concatenate((data,data,data),-1)
    else:
        data = N.concatenate((data,data,data),-1)

    lon = lon - 360
    # clip grid - this allows to clip map anywhere between -360 and 360
    ind1 = N.searchsorted(lon[0,:],-360)
    ind2 = N.searchsorted(lon[0,:],360)
    lon  = lon[:,ind1-1:ind2+1]
    lat  = lat[:,ind1-1:ind2+1]
    data = data[...,ind1-1:ind2+1]

    return lon, lat, data

def _map_common(draw_method_name, arg_func, mode, cube, data, *args, **kwargs):
    """
    Draw the given cube on a map using its points or bounds.

    "Mode" parameter will switch functionality between POINT or BOUND plotting.

    """
    # get the 2d x and 2d y from the CS
    if mode == iris.coords.POINT_MODE:
        x, y = cartography.get_xy_grids(cube)
    else:
        try:
            x, y = cartography.get_xy_contiguous_bounded_grids(cube)
        # Exception translation.
        except iris.exceptions.CoordinateMultiDimError:
            raise ValueError("Could not get XY grid from bounds. "
                             "X or Y coordinate not 1D.")
        except ValueError:
            raise ValueError("Could not get XY grid from bounds. "
                             "X or Y coordinate doesn't have 2 bounds "
                             "per point.")

    # take a copy of the data so that we can make modifications to it
    data = data.copy()

    # If we are global, then append the first column of data the array to the
    # last (and add 360 degrees) NOTE: if it is found that this block of code
    # is useful in anywhere other than this plotting routine, it may be better
    # placed in the CS.
    x_coord = cube.coord(axis="X")
    if getattr(x_coord, 'circular', False):
        _, direction = iris.util.monotonic(x_coord.points,
                                           return_direction=True)
        y = np.append(y, y[:, 0:1], axis=1)
        x = np.append(x, x[:, 0:1] + 360 * direction, axis=1)
        data = ma.concatenate([data, data[:, 0:1]], axis=1)

    # Replace non-cartopy subplot/axes with a cartopy alternative.
    cs = cube.coord_system('CoordSystem')
    if cs:
        cartopy_proj = cs.as_cartopy_projection()
    else:
        cartopy_proj = cartopy.crs.PlateCarree()
    ax = _get_cartopy_axes(cartopy_proj)

    draw_method = getattr(ax, draw_method_name)

    # Set the "from transform" keyword.
    # NB. While cartopy doesn't support spherical contours, just use the
    # projection as the source CRS.
    assert 'transform' not in kwargs, 'Transform keyword is not allowed.'
    kwargs['transform'] = cartopy_proj

    if arg_func is not None:
        new_args, kwargs = arg_func(x, y, data, *args, **kwargs)
    else:
        new_args = (x, y, data) + args

    # Draw the contour lines/filled contours.
    return draw_method(*new_args, **kwargs)

def _map_common(draw_method_name, arg_func, mode, cube, plot_defn,
                *args, **kwargs):
    """
    Draw the given cube on a map using its points or bounds.

    "Mode" parameter will switch functionality between POINT or BOUND plotting.


    """
    # Generate 2d x and 2d y grids.
    y_coord, x_coord = plot_defn.coords
    if mode == iris.coords.POINT_MODE:
        if x_coord.ndim == y_coord.ndim == 1:
            x, y = np.meshgrid(x_coord.points, y_coord.points)
        elif x_coord.ndim == y_coord.ndim == 2:
            x = x_coord.points
            y = y_coord.points
        else:
            raise ValueError("Expected 1D or 2D XY coords")
    else:
        try:
            x, y = np.meshgrid(x_coord.contiguous_bounds(),
                               y_coord.contiguous_bounds())
        # Exception translation.
        except iris.exceptions.CoordinateMultiDimError:
            raise ValueError("Could not get XY grid from bounds. "
                             "X or Y coordinate not 1D.")
        except ValueError:
            raise ValueError("Could not get XY grid from bounds. "
                             "X or Y coordinate doesn't have 2 bounds "
                             "per point.")

    # Obtain the data array.
    data = cube.data
    if plot_defn.transpose:
        data = data.T

    # If we are global, then append the first column of data the array to the
    # last (and add 360 degrees) NOTE: if it is found that this block of code
    # is useful in anywhere other than this plotting routine, it may be better
    # placed in the CS.
    if getattr(x_coord, 'circular', False):
        _, direction = iris.util.monotonic(x_coord.points,
                                           return_direction=True)
        y = np.append(y, y[:, 0:1], axis=1)
        x = np.append(x, x[:, 0:1] + 360 * direction, axis=1)
        data = ma.concatenate([data, data[:, 0:1]], axis=1)

    # Replace non-cartopy subplot/axes with a cartopy alternative and set the
    # transform keyword.
    kwargs = _ensure_cartopy_axes_and_determine_kwargs(x_coord, y_coord,
                                                       kwargs)

    if arg_func is not None:
        new_args, kwargs = arg_func(x, y, data, *args, **kwargs)
    else:
        new_args = (x, y, data) + args

    # Draw the contour lines/filled contours.
    return getattr(plt, draw_method_name)(*new_args, **kwargs)

 def extend_interp(datafield):
   # add masked values at southernmost end
   southernlimitmask = ma.masked_all(len(self.olon))
   olat_ext          = np.append(-82.1,self.olat)
   dfield_ext = ma.concatenate([ma.column_stack(southernlimitmask), datafield], 0)
   # f = interp2d(self.olon, olat_ext, dfield_ext)
   # return f(self.pismlon, self.pismlat)
   return interp(dfield_ext, self.olon, olat_ext, self.pismlon, self.pismlat)

    def recache(self, always=False):
        if always or self._invalidx:
            xconv = self.convert_xunits(self._xorig)
            if ma.isMaskedArray(self._xorig):
                x = ma.asarray(xconv, np.float_)
            else:
                x = np.asarray(xconv, np.float_)
            x = x.ravel()
        else:
            x = self._x
        if always or self._invalidy:
            yconv = self.convert_yunits(self._yorig)
            if ma.isMaskedArray(self._yorig):
                y = ma.asarray(yconv, np.float_)
            else:
                y = np.asarray(yconv, np.float_)
            y = y.ravel()
        else:
            y = self._y

        if len(x) == 1 and len(y) > 1:
            x = x * np.ones(y.shape, np.float_)
        if len(y) == 1 and len(x) > 1:
            y = y * np.ones(x.shape, np.float_)

        if len(x) != len(y):
            raise RuntimeError("xdata and ydata must be the same length")

        x = x.reshape((len(x), 1))
        y = y.reshape((len(y), 1))

        if ma.isMaskedArray(x) or ma.isMaskedArray(y):
            self._xy = ma.concatenate((x, y), 1)
        else:
            self._xy = np.concatenate((x, y), 1)
        self._x = self._xy[:, 0]  # just a view
        self._y = self._xy[:, 1]  # just a view

        self._subslice = False
        if (
            self.axes
            and len(x) > 100
            and self._is_sorted(x)
            and self.axes.name == "rectilinear"
            and self.axes.get_xscale() == "linear"
            and self._markevery is None
        ):
            self._subslice = True
        if hasattr(self, "_path"):
            interpolation_steps = self._path._interpolation_steps
        else:
            interpolation_steps = 1
        self._path = Path(self._xy, None, interpolation_steps)
        self._transformed_path = None
        self._invalidx = False
        self._invalidy = False

def test_flags2bin(n=100):
    flag = ma.concatenate([np.random.randint(0,5,n),
        ma.masked_all(2, dtype='int8')])

    binflags = flags2bin(flag)

    assert type(binflags) == ma.MaskedArray
    assert binflags.dtype == 'bool'
    assert binflags.shape == (n+2,)
    assert binflags.mask[flag.mask].all(), \
            "All masked flags records should be also masked at binflags"

    def test_testCopySize(self):
        # Tests of some subtle points of copying and sizing.
        n = [0, 0, 1, 0, 0]
        m = make_mask(n)
        m2 = make_mask(m)
        assert_(m is m2)
        m3 = make_mask(m, copy=1)
        assert_(m is not m3)

        x1 = np.arange(5)
        y1 = array(x1, mask=m)
        assert_(y1._data is not x1)
        assert_(allequal(x1, y1._data))
        assert_(y1._mask is m)

        y1a = array(y1, copy=0)
        # For copy=False, one might expect that the array would just
        # passed on, i.e., that it would be "is" instead of "==".
        # See gh-4043 for discussion.
        assert_(y1a._mask.__array_interface__ ==
                y1._mask.__array_interface__)

        y2 = array(x1, mask=m3, copy=0)
        assert_(y2._mask is m3)
        assert_(y2[2] is masked)
        y2[2] = 9
        assert_(y2[2] is not masked)
        assert_(y2._mask is m3)
        assert_(allequal(y2.mask, 0))

        y2a = array(x1, mask=m, copy=1)
        assert_(y2a._mask is not m)
        assert_(y2a[2] is masked)
        y2a[2] = 9
        assert_(y2a[2] is not masked)
        assert_(y2a._mask is not m)
        assert_(allequal(y2a.mask, 0))

        y3 = array(x1 * 1.0, mask=m)
        assert_(filled(y3).dtype is (x1 * 1.0).dtype)

        x4 = arange(4)
        x4[2] = masked
        y4 = resize(x4, (8,))
        assert_(eq(concatenate([x4, x4]), y4))
        assert_(eq(getmask(y4), [0, 0, 1, 0, 0, 0, 1, 0]))
        y5 = repeat(x4, (2, 2, 2, 2), axis=0)
        assert_(eq(y5, [0, 0, 1, 1, 2, 2, 3, 3]))
        y6 = repeat(x4, 2, axis=0)
        assert_(eq(y5, y6))

 def append(self, na):
     """
 Append tableDict na to self. 
 """
     # for key in na.cols:
     # if key == 'keys': continue
     # if key in self.cols:
     # newDataVec = self.__getitem__(key)
     # for ele in na.data[key]:
     ##newDataVec.append(ele)
     # np.append(newDataVec, ele)
     # self.__setitem__(key, newDataVec)
     for col in self.cols:
         self.data[col] = ma.concatenate((self.data[col][:], na.data[col][:]))
     self.nRows = self.nRows + na.nRows

def densitystep(S, T, P):
    """
    """
    assert S.shape == T.shape
    assert S.shape == P.shape
    try:
        import gsw
        rho0 = gsw.pot_rho_t_exact(S, T, P, 0)
        assert S.ndim == 1, "Not able to densitystep an array ndim > 1"
        ds = ma.concatenate([ma.masked_all(1),
                np.sign(np.diff(P))*np.diff(rho0)])
        return ma.fix_invalid(ds)

    except ImportError:
        print("Package gsw is required and is not available.")

def addcyclic(data):
	"""
	Adds cyclic points to an array in rightmost dimension.
	data = input 2D array.
	"""
	if data.ndim != 2:
		print('ERROR: Input array is not two-dimensional')
		return

	if MA.isMA(data):
		newdata = MA.concatenate((data,data[:,0,N.newaxis]),axis=-1)
	else:
		newdata = N.concatenate((data,data[:,0,N.newaxis]),axis=-1)

	return newdata

def fix_pop_grid_test(tlon,tlat,data):
    """
    Pad coordinates and data on CCSM/POP gx3v5 grid so it can be
    plotted with matplotlib/basemap
    tlon,tlat,data must be 2D arrays
    """
    # make lon monotonic and pad coordinate and data arrays along lon axis
    tlon = N.where(N.greater_equal(tlon,min(tlon[:,0])),tlon-360,tlon)
    lon  = N.concatenate((tlon,tlon+360),1)
    lat  = N.concatenate((tlat,tlat),1)
    if hasattr(data,'mask'):
        data = MA.concatenate((data,data),1)
    else:
        data = N.concatenate((data,data),1)

    return lon, lat, data

def test_i2b_flags(n=100):
    flag = ma.concatenate([np.random.randint(0,5,n),
        ma.masked_all(2, dtype='int8')])

    binflags = i2b_flags(flag)

    assert type(binflags) == ma.MaskedArray
    assert binflags.dtype == 'bool'
    assert binflags.shape == (n+2,)
    assert binflags.mask[flag.mask].all(), \
            "All masked flags records should be also masked at binflags"

    # FIXME: Improve this. Include cases with dict as input.
    #   Check the output in case of dict input and differente combinations
    #     of True/False, masked or not.
    assert (i2b_flags([1, 2, 3, 4]) ==
            ma.array([True, True, False, False])).all()

    def _build_data(self):
        """
        Generate the data payload for the new concatenated cube.

        Returns:
            The concatenated :class:`iris.cube.Cube` data payload.

        """
        skeletons = self._skeletons
        data = [skeleton.data for skeleton in skeletons]

        if self._data_is_masked:
            data = ma.concatenate(tuple(data), axis=self.axis)
        else:
            data = np.concatenate(tuple(data), axis=self.axis)

        return data

    def test_testCopySize(self):
        # Tests of some subtle points of copying and sizing.
        n = [0, 0, 1, 0, 0]
        m = make_mask(n)
        m2 = make_mask(m)
        assert_(m is m2)
        m3 = make_mask(m, copy=1)
        assert_(m is not m3)

        x1 = np.arange(5)
        y1 = array(x1, mask=m)
        assert_(y1._data is not x1)
        assert_(allequal(x1, y1._data))
        assert_(y1.mask is m)

        y1a = array(y1, copy=0)
        assert_(y1a.mask is y1.mask)

        y2 = array(x1, mask=m3, copy=0)
        assert_(y2.mask is m3)
        assert_(y2[2] is masked)
        y2[2] = 9
        assert_(y2[2] is not masked)
        assert_(y2.mask is m3)
        assert_(allequal(y2.mask, 0))

        y2a = array(x1, mask=m, copy=1)
        assert_(y2a.mask is not m)
        assert_(y2a[2] is masked)
        y2a[2] = 9
        assert_(y2a[2] is not masked)
        assert_(y2a.mask is not m)
        assert_(allequal(y2a.mask, 0))

        y3 = array(x1 * 1.0, mask=m)
        assert_(filled(y3).dtype is (x1 * 1.0).dtype)

        x4 = arange(4)
        x4[2] = masked
        y4 = resize(x4, (8,))
        assert_(eq(concatenate([x4, x4]), y4))
        assert_(eq(getmask(y4), [0, 0, 1, 0, 0, 0, 1, 0]))
        y5 = repeat(x4, (2, 2, 2, 2), axis=0)
        assert_(eq(y5, [0, 0, 1, 1, 2, 2, 3, 3]))
        y6 = repeat(x4, 2, axis=0)
        assert_(eq(y5, y6))

 def test_testAddSumProd(self):
     # Test add, sum, product.
     (x, y, a10, m1, m2, xm, ym, z, zm, xf, s) = self.d
     self.assertTrue(eq(np.add.reduce(x), add.reduce(x)))
     self.assertTrue(eq(np.add.accumulate(x), add.accumulate(x)))
     self.assertTrue(eq(4, sum(array(4), axis=0)))
     self.assertTrue(eq(4, sum(array(4), axis=0)))
     self.assertTrue(eq(np.sum(x, axis=0), sum(x, axis=0)))
     self.assertTrue(eq(np.sum(filled(xm, 0), axis=0), sum(xm, axis=0)))
     self.assertTrue(eq(np.sum(x, 0), sum(x, 0)))
     self.assertTrue(eq(np.product(x, axis=0), product(x, axis=0)))
     self.assertTrue(eq(np.product(x, 0), product(x, 0)))
     self.assertTrue(eq(np.product(filled(xm, 1), axis=0), product(xm, axis=0)))
     if len(s) > 1:
         self.assertTrue(eq(np.concatenate((x, y), 1), concatenate((xm, ym), 1)))
         self.assertTrue(eq(np.add.reduce(x, 1), add.reduce(x, 1)))
         self.assertTrue(eq(np.sum(x, 1), sum(x, 1)))
         self.assertTrue(eq(np.product(x, 1), product(x, 1)))