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_).filled(np.nan)
            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_).filled(np.nan)
            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')

        self._xy = np.empty((len(x), 2), dtype=np.float_)
        self._xy[:, 0] = x
        self._xy[:, 1] = y

        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) > 1000 and self._is_sorted(x) and
                self.axes.name == 'rectilinear' and
                self.axes.get_xscale() == 'linear' and
                self._markevery is None and
                self.get_clip_on() is True):
            self._subslice = True
            nanmask = np.isnan(x)
            if nanmask.any():
                self._x_filled = self._x.copy()
                indices = np.arange(len(x))
                self._x_filled[nanmask] = np.interp(indices[nanmask],
                        indices[~nanmask], self._x[~nanmask])
            else:
                self._x_filled = self._x

        if self._path is not None:
            interpolation_steps = self._path._interpolation_steps
        else:
            interpolation_steps = 1
        xy = STEP_LOOKUP_MAP[self._drawstyle](*self._xy.T)
        self._path = Path(np.asarray(xy).T, None, interpolation_steps)
        self._transformed_path = None
        self._invalidx = False
        self._invalidy = False

def to_ham6(img, palette, background=None, out=None):
    _debug_array(img)

    if background is None:
        background = ma.masked
    elif isinstance(background, numbers.Integral):
        background = palette[background]

    if not ma.is_masked(background) and ma.isMaskedArray(img):
        img = img.filled(background)

    if ma.isMaskedArray(img):
        ham6 = ma.empty(img.shape[:2], dtype=np.uint8)
    else:
        ham6 = np.empty(img.shape[:2], dtype=np.uint8)

    for y in range(img.shape[0]):
        c = background
        for x in range(img.shape[1]):
            i, c = ham6_nearest(img[y, x], palette, c)
            ham6[y, x] = i
            if out is not None:
                out[y, x] = c

    _debug_array(ham6)
    return ham6

def generic_interp_pres(p, pres, field):
    '''
    Generic interpolation routine

    Parameters
    ----------
    p : number, numpy array
        Pressure (hPa) of the level for which the field variable is desired
    pres : numpy array
        The array of pressure
    field : numpy array
        The variable which is being interpolated
    log : bool
        Flag to determine whether the 'field' variable is in log10 space

    Returns
    -------
    Value of the 'field' variable at the given pressure

    '''
    if ma.isMaskedArray(pres):
        not_masked1 = ~pres.mask
    else:
        not_masked1 = np.ones(pres.shape, dtype=bool)
        not_masked1[:] = True
    if ma.isMaskedArray(field):
        not_masked2 = ~field.mask
    else:
        not_masked2 = np.ones(field.shape, dtype=bool)
        not_masked2[:] = True
    not_masked = not_masked1 * not_masked2
    return np.interp(p, pres[not_masked], field[not_masked], left=ma.masked,
                     right=ma.masked)

    def test_peak_with_mask(self):
        # Single value in column masked.
        latitude = iris.coords.DimCoord(np.arange(0, 5, 1),
                                        standard_name='latitude',
                                        units='degrees')
        cube = iris.cube.Cube(ma.array([1, 4, 2, 3, 2], dtype=np.float32),
                              standard_name='air_temperature',
                              units='kelvin')
        cube.add_dim_coord(latitude, 0)

        cube.data[3] = ma.masked

        collapsed_cube = cube.collapsed('latitude', iris.analysis.PEAK)
        self.assertArrayAlmostEqual(collapsed_cube.data,
                                    np.array([4.024977], dtype=np.float32))
        self.assertTrue(ma.isMaskedArray(collapsed_cube.data))
        self.assertEqual(collapsed_cube.data.shape, (1,))

        # Whole column masked.
        cube.data[:] = ma.masked

        collapsed_cube = cube.collapsed('latitude', iris.analysis.PEAK)
        masked_array = ma.array(ma.masked)
        self.assertTrue(ma.allequal(collapsed_cube.data, masked_array))
        self.assertTrue(ma.isMaskedArray(collapsed_cube.data))
        self.assertEqual(collapsed_cube.data.shape, (1,))

    def test_peak_with_nan_and_mask(self):
        # Single nan in column with single value masked.
        latitude = iris.coords.DimCoord(np.arange(0, 5, 1),
                                        standard_name='latitude',
                                        units='degrees')
        cube = iris.cube.Cube(ma.array([1, 4, 2, 3, 1], dtype=np.float32),
                              standard_name='air_temperature',
                              units='kelvin')
        cube.add_dim_coord(latitude, 0)

        cube.data[3] = np.nan
        cube.data[4] = ma.masked

        collapsed_cube = cube.collapsed('latitude', iris.analysis.PEAK)
        self.assertArrayAlmostEqual(collapsed_cube.data,
                                    np.array([4.024977], dtype=np.float32))
        self.assertTrue(ma.isMaskedArray(collapsed_cube.data))
        self.assertEqual(collapsed_cube.data.shape, (1,))

        # Only nans in column where values not masked.
        cube.data[0:3] = np.nan

        collapsed_cube = cube.collapsed('latitude', iris.analysis.PEAK)
        self.assertTrue(np.isnan(collapsed_cube.data).all())
        self.assertTrue(ma.isMaskedArray(collapsed_cube.data))
        self.assertEqual(collapsed_cube.data.shape, (1,))

    def test_rotated_to_osgb(self):
        # Rotated Pole data with large extent.
        x = np.linspace(311.9, 391.1, 10)
        y = np.linspace(-23.6, 24.8, 8)
        u, v = uv_cubes(x, y)
        ut, vt = rotate_winds(u, v, iris.coord_systems.OSGB())

        # Ensure cells with discrepancies in magnitude are masked.
        self.assertTrue(ma.isMaskedArray(ut.data))
        self.assertTrue(ma.isMaskedArray(vt.data))

        # Snapshot of mask with fixed tolerance of atol=2e-3
        expected_mask = np.array([[1, 1, 1, 0, 0, 0, 0, 0, 0, 1],
                                  [1, 1, 1, 0, 0, 0, 0, 0, 0, 1],
                                  [1, 1, 1, 1, 0, 0, 0, 0, 1, 1],
                                  [1, 1, 1, 1, 0, 0, 0, 0, 1, 1],
                                  [1, 1, 1, 1, 0, 0, 0, 0, 1, 1],
                                  [1, 1, 1, 1, 1, 0, 0, 1, 1, 1],
                                  [1, 1, 1, 1, 1, 0, 0, 1, 1, 1],
                                  [1, 1, 1, 1, 1, 0, 0, 1, 1, 1]], np.bool)
        self.assertArrayEqual(expected_mask, ut.data.mask)
        self.assertArrayEqual(expected_mask, vt.data.mask)

        # Check unmasked values have sufficiently small error in mag.
        expected_mag = np.sqrt(u.data**2 + v.data**2)
        # Use underlying data to ignore mask in calculation.
        res_mag = np.sqrt(ut.data.data**2 + vt.data.data**2)
        # Calculate percentage error (note there are no zero magnitudes
        # so we can divide safely).
        anom = 100.0 * np.abs(res_mag - expected_mag) / expected_mag
        self.assertTrue(anom[~ut.data.mask].max() < 0.1)

def generic_interp_hght(h, hght, field, log=False):
    '''
    Generic interpolation routine

    Parameters
    ----------
    h : number, numpy array
        Height (m) of the level for which pressure is desired
    hght : numpy array
        The array of heights
    field : numpy array
        The variable which is being interpolated
    log : bool
        Flag to determine whether the 'field' variable is in log10 space

    Returns
    -------
    Value of the 'field' variable at the given height

    '''
    if ma.isMaskedArray(hght):
        not_masked1 = ~hght.mask
    else:
        not_masked1 = np.ones(hght.shape)
    if ma.isMaskedArray(field):
        not_masked2 = ~field.mask
    else:
        not_masked2 = np.ones(field.shape)
    not_masked = not_masked1 * not_masked2
    if log:
        return 10**np.interp(h, hght[not_masked], field[not_masked],
                             left=ma.masked, right=ma.masked)
    else:
        return np.interp(h, hght[not_masked], field[not_masked],
                         left=ma.masked, right=ma.masked)

    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 __init__(self, cpfile, smapfile, areafile, crop, varname):
        with nc(cpfile) as f:
            self.year    = f.variables['year'][:]
            self.week    = f.variables['week'][:]
            self.county  = f.variables['county'][:]
            self.day     = f.variables['day'][:]
            self.rawdata = f.variables[varname][:]

            varatt = f.variables['var'].ncattrs()
            if 'units' in varatt and f.variables['var'].units == 'mapping':
                self.var = array(f.variables['var'].long_name.split(', '))
                self.varmap = self.varmap_str[crop]
            else:          
                self.var = f.variables['var'][:]
                self.varmap = self.varmap_num[crop]

        self.crop = crop

        nyears, nweeks, ncounties, nvars = len(self.year), len(self.week), len(self.county), len(self.vars)

        self.data = masked_array(zeros((nyears, nweeks, ncounties, nvars)), mask = ones((nyears, nweeks, ncounties, nvars)))
        for i in range(nvars):
            vmap = self.varmap[i]

            if isinstance(vmap, list):
                for j in range(ncounties):
                    for k in range(len(vmap)):
                        if vmap[k] in self.var: # variable in list
                            varidx = where(self.var == vmap[k])[0][0]
                            data   = self.rawdata[:, :, varidx, j]
                            if not isMaskedArray(data) or not data.mask.all():
                                self.data[:, :, j, i] = data
                                break
            elif vmap != '':
                if vmap in self.var:
                    varidx = where(self.var == vmap)[0][0]
                    self.data[:, :, :, i] = self.rawdata[:, :, varidx, :]
            else: # no data
                continue

            # discard counties with insufficient data
            for j in range(ncounties):
                for k in range(nyears):
                    data = self.data[k, :, j, i]

                    if isMaskedArray(data):
                        data = data[~data.mask]

                    if data.size and data[-1] - data[0] < 40:
                        self.data[k, :, j, i].mask = True # mask

        # aggregate to state level
        aggregator = StateAggregator(smapfile, areafile)
        self.sdata = aggregator.aggregate(self.data, self.county)
        self.state = aggregator.states

 def _check_fill_value(self, result, fill0='none', fill1='none'):
     expected_fill_value = self._expected_fill_value(fill0, fill1)
     if expected_fill_value is None:
         data = result.data
         if ma.isMaskedArray(data):
             np_fill_value = ma.masked_array(0,
                                             dtype=result.dtype).fill_value
             self.assertEqual(data.fill_value, np_fill_value)
     else:
         data = result.data
         if ma.isMaskedArray(data):
             self.assertEqual(data.fill_value, expected_fill_value)

def generic_interp_hght(h, hght, field, log=False):
    '''
    Generic interpolation routine

    Parameters
    ----------
    h : number, numpy array
        Height (m) of the level for which pressure is desired
    hght : numpy array
        The array of heights
    field : numpy array
        The variable which is being interpolated
    log : bool
        Flag to determine whether the 'field' variable is in log10 space

    Returns
    -------
    Value of the 'field' variable at the given height

    '''
    if ma.isMaskedArray(hght):
        not_masked1 = ~hght.mask
    else:
        not_masked1 = np.ones(hght.shape)
    if ma.isMaskedArray(field):
        not_masked2 = ~field.mask
    else:
        not_masked2 = np.ones(field.shape)
    not_masked = not_masked1 * not_masked2

    field_intrp = np.interp(h, hght[not_masked], field[not_masked],
                         left=ma.masked, right=ma.masked)
 
    if hasattr(h, 'shape') and h.shape == tuple():
        h = h[()]

    if type(h) != type(ma.masked):
        # Bug fix for Numpy v1.10: returns nan on the boundary.
        field_intrp = np.where(np.isclose(h, hght[not_masked][0]), field[not_masked][0], field_intrp)
        field_intrp = np.where(np.isclose(h, hght[not_masked][-1]), field[not_masked][-1], field_intrp)

    # Another bug fix: np.interp() returns masked values as nan. We want ma.masked, dangit!
    field_intrp = ma.where(np.isnan(field_intrp), ma.masked, field_intrp)

    # ma.where() returns a 0-d array when the arguments are floats, which confuses subsequent code.
    if hasattr(field_intrp, 'shape') and field_intrp.shape == tuple():
        field_intrp = field_intrp[()]

    if log:
        return 10 ** field_intrp
    else:
        return field_intrp

 def test_testBasic1d(self):
     # Test of basic array creation and properties in 1 dimension.
     (x, y, a10, m1, m2, xm, ym, z, zm, xf, s) = self.d
     assert_(not isMaskedArray(x))
     assert_(isMaskedArray(xm))
     assert_equal(shape(xm), s)
     assert_equal(xm.shape, s)
     assert_equal(xm.dtype, x.dtype)
     assert_equal(xm.size, reduce(lambda x, y:x * y, s))
     assert_equal(count(xm), len(m1) - reduce(lambda x, y:x + y, m1))
     assert_(eq(xm, xf))
     assert_(eq(filled(xm, 1.e20), xf))
     assert_(eq(x, xm))

def generic_interp_pres(p, pres, field):
    '''
    Generic interpolation routine

    Parameters
    ----------
    p : number, numpy array
        Pressure (hPa) of the level for which the field variable is desired
    pres : numpy array
        The array of pressure
    field : numpy array
        The variable which is being interpolated
    log : bool
        Flag to determine whether the 'field' variable is in log10 space

    Returns
    -------
    Value of the 'field' variable at the given pressure

    '''
    if ma.isMaskedArray(pres):
        not_masked1 = ~pres.mask
    else:
        not_masked1 = np.ones(pres.shape, dtype=bool)
        not_masked1[:] = True
    if ma.isMaskedArray(field):
        not_masked2 = ~field.mask
    else:
        not_masked2 = np.ones(field.shape, dtype=bool)
        not_masked2[:] = True
    not_masked = not_masked1 * not_masked2

    field_intrp = np.interp(p, pres[not_masked], field[not_masked], left=ma.masked,
                 right=ma.masked)

    if hasattr(p, 'shape') and p.shape == tuple():
        p = p[()]

    if type(p) != type(ma.masked):
        # Bug fix for Numpy v1.10: returns nan on the boundary.
        field_intrp = ma.where(np.isclose(p, pres[not_masked][0]), field[not_masked][0], field_intrp)
        field_intrp = ma.where(np.isclose(p, pres[not_masked][-1]), field[not_masked][-1], field_intrp)

    # Another bug fix: np.interp() returns masked values as nan. We want ma.masked, dangit!
    field_intrp = ma.where(np.isnan(field_intrp), ma.masked, field_intrp)

    # ma.where() returns a 0-d array when the arguments are floats, which confuses subsequent code.
    if hasattr(field_intrp, 'shape') and field_intrp.shape == tuple():
        field_intrp = field_intrp[()]

    return field_intrp

def _weighted_mean_with_mdtol(data, weights, axis=None, mdtol=0):
    """
    Return the weighted mean of an array over the specified axis
    using the provided weights (if any) and a permitted fraction of
    masked data.

    Args:

    * data (array-like):
        Data to be averaged.

    * weights (array-like):
        An array of the same shape as the data that specifies the contribution
        of each corresponding data element to the calculated mean.

    Kwargs:

    * axis (int or tuple of ints):
        Axis along which the mean is computed. The default is to compute
        the mean of the flattened array.

    * mdtol (float):
        Tolerance of missing data. The value returned in each element of the
        returned array will be masked if the fraction of masked data exceeds
        mdtol. This fraction is weighted by the `weights` array if one is
        provided. mdtol=0 means no missing data is tolerated
        while mdtol=1 will mean the resulting element will be masked if and
        only if all the contributing elements of data are masked.
        Defaults to 0.

    Returns:
        Numpy array (possibly masked) or scalar.

    """
    res = ma.average(data, weights=weights, axis=axis)
    if ma.isMaskedArray(data) and mdtol < 1:
        weights_total = weights.sum(axis=axis)
        masked_weights = weights.copy()
        masked_weights[~ma.getmaskarray(data)] = 0
        masked_weights_total = masked_weights.sum(axis=axis)
        frac_masked = np.true_divide(masked_weights_total, weights_total)
        mask_pt = frac_masked > mdtol
        if np.any(mask_pt):
            if np.isscalar(res):
                res = ma.masked
            elif ma.isMaskedArray(res):
                res.mask |= mask_pt
            else:
                res = ma.masked_array(res, mask=mask_pt)
    return res

    def __loadwts(self, var, vals):
        if self.wfile is None:
            return ones(len(vals))
        else:
            vars  = var.split('/')
            nvars = len(vars)

            v = [0] * nvars
            w = [0] * nvars
            with nc(self.wfile) as f:
                for i in range(nvars):
                    if f.variables[vars[i]].units == 'mapping':
                        v[i] = array(f.variables[vars[i]].long_name.split(', '))
                    else:
                        v[i] = f.variables[vars[i]][:]
                    w[i] = f.variables['weights_' + vars[i]][:]

            nvals = len(vals)
            wts = masked_array(zeros(nvals), mask = ones(nvals))
            for i in range(nvals):
                svals = vals[i].split('/')
                for j in range(nvars):
                    if svals[j].isdigit():
                        svals[j] = double(svals[j]) # convert to number
                    idx = where(v[j] == svals[j])[0]
                    if idx.size:
                        if isMaskedArray(wts[i]):
                            wts[i] = w[j][idx[0]]
                        else:
                            wts[i] *= w[j][idx[0]]

            return wts

def concatenate(cubes, order=None):
    """
    Explicitly force the contiguous major order of cube data
    alignment to ensure consistent CML crc32 checksums.

    Defaults to contiguous 'C' row-major order.

    """
    if order is None:
        order = 'C'

    cubelist = iris.cube.CubeList(cubes)
    result = cubelist.concatenate()

    for cube in result:
        # Setting the cube.data clears the cube.dtype and cube.fill_value.
        # We want to maintain the fill_value for testing purposes, even
        # though we have lost the lazy data in order to pin down the array
        # data order to overcome testing on different architectures.
        fill_value = cube.fill_value
        if ma.isMaskedArray(cube.data):
            data = np.array(cube.data.data, copy=True, order=order)
            mask = np.array(cube.data.mask, copy=True, order=order)
            cube.data = ma.array(data, mask=mask)
        else:
            cube.data = np.array(cube.data, copy=True, order=order)
        cube.fill_value = fill_value

    return result

    def identify(cls, variables, ignore=None, target=None, warn=True, monotonic=False):
        result = {}
        ignore, target = cls._identify_common(variables, ignore, target)

        # Identify all CF coordinate variables.
        for nc_var_name, nc_var in target.iteritems():
            if nc_var_name in ignore:
                continue
            # String variables can't be coordinates
            if np.issubdtype(nc_var.dtype, np.str):
                continue
            # Restrict to one-dimensional with name as dimension OR zero-dimensional scalar
            if not ((nc_var.ndim == 1 and nc_var_name in nc_var.dimensions) or (nc_var.ndim == 0)):
                continue
            # Restrict to monotonic?
            if monotonic:
                data = nc_var[:]
                # Gracefully fill a masked coordinate.
                if ma.isMaskedArray(data):
                    data = ma.filled(data)
                if nc_var.shape == () or nc_var.shape == (1,) or iris.util.monotonic(data):
                    result[nc_var_name] = CFCoordinateVariable(nc_var_name, nc_var)
            else:
                result[nc_var_name] = CFCoordinateVariable(nc_var_name, nc_var)

        return result

def _condition_arg(value):
    """
    Validate value is acceptable for conversion purposes.

    Will convert into an array if not a scalar, and can be converted
    into an array

    Parameters
    ----------
    value: int or float value, or sequence of such values

    Returns
    -------
    Scalar value or numpy array

    Raises
    ------
    ValueError
        If value is not as expected
    """
    if isinstance(value, (float, int, long)):
        return value
    else:
        try:
            avalue = np.array(value)
            dt = str(avalue.dtype)
            if not (dt.startswith('int') or dt.startswith('float')):
                raise ValueError("Must be convertable to int or float array")
            if ma.isMaskedArray(value):
                return value
            return avalue
        except ValueError:
            raise ValueError(
                "Value not scalar compatible or convertable into a float or "
                "integer array")

 def test_scalar_no_overlap(self):
     # Slice src so result collapses to a scalar.
     src_cube = self.src_cube[:, 1, :]
     # Regrid to a single cell with no overlap with masked src cells.
     grid_cube = self.grid_cube[2, 1, 3]
     res = regrid(src_cube, grid_cube, mdtol=0.8)
     self.assertFalse(ma.isMaskedArray(res.data))

def _get_scale_factor_add_offset(cube, datatype):
    """Utility function used by netCDF data packing tests."""
    if isinstance(datatype, dict):
        dt = np.dtype(datatype['dtype'])
    else:
        dt = np.dtype(datatype)
    cmax = cube.data.max()
    cmin = cube.data.min()
    n = dt.itemsize * 8
    if ma.isMaskedArray(cube.data):
        masked = True
    else:
        masked = False
    if masked:
        scale_factor = (cmax - cmin)/(2**n-2)
    else:
        scale_factor = (cmax - cmin)/(2**n-1)
    if dt.kind == 'u':
        add_offset = cmin
    elif dt.kind == 'i':
        if masked:
            add_offset = (cmax + cmin)/2
        else:
            add_offset = cmin + 2**(n-1)*scale_factor
    return (scale_factor, add_offset)