def make_gene_map_1(self):
        count = 0
        self.iterator = itertools.product(range(self.time_len), range(self.lat_len), range(self.lon_len))
        for x_valid in self.iterator:
            binary_string = bin(count)[2:]
            while len(binary_string) < self.string_length:
                binary_string = "0" + binary_string
                # self.gene_map[binary_string] = {}
                # self.gene_map[binary_string]['coordinate'] = tuple(x_valid)
            if ma.getdata(self.array[x_valid]) < 100:  # chooses unmasked points?
                self.gene_map[binary_string] = {}
                self.gene_map[binary_string]["coordinate"] = tuple(x_valid)
                self.gene_map[binary_string]["value"] = self.array[x_valid]
                # count += 1
                print count, binary_string, self.gene_map[binary_string]["value"]  # debug
            else:
                # self.gene_map[binary_string]['value'] = 1E06
                pass
            count += 1
            # print count, binary_string, self.gene_map[binary_string]['value'] # debug
        self.last_valid_binary_string = binary_string
        not_valid_first = eval("0b" + binary_string) + 1
        not_valid_last = eval("0b" + "1" * self.string_length)
        print not_valid_last
        print x_valid
        print ma.isMA(self.array)

        for x_not_valid in range(not_valid_first, not_valid_last + 1):
            binary_string = bin(x_not_valid)[2:]
            self.gene_map[binary_string] = {}
            self.gene_map[binary_string]["coordinate"] = (999, 999, 999)
            self.gene_map[binary_string]["value"] = 1e06
            print x_not_valid, binary_string, self.gene_map[binary_string]["value"]

def shiftgrid(lon0, datain, lonsin, start= True, cyclic=360.0):
    """
    Purpose::
        Shift global lat/lon grid east or west. This function is taken directly
        from the (unreleased) basemap 1.0.7 source code as version 1.0.6 does not
        currently support arrays with more than two dimensions.
        https://github.com/matplotlib/basemap
        
    Input::
        lon0 - starting longitude for shifted grid (ending longitude if start=False). 
               lon0 must be on input grid (within the range of lonsin).
        datain - original data with longitude the right-most dimension.
        lonsin - original longitudes.
        start  - if True, lon0 represents the starting longitude of the new grid. 
                 if False, lon0 is the ending longitude. Default True.
        cyclic - width of periodic domain (default 360)

    Output:: 
        dataout - data on shifted grid
        lonsout - lons on shifted grid
    """
    if np.fabs(lonsin[-1]-lonsin[0]-cyclic) > 1.e-4:
        # Use all data instead of raise ValueError, 'cyclic point not included'
        start_idx = 0
    else:
        # If cyclic, remove the duplicate point
        start_idx = 1
    if lon0 < lonsin[0] or lon0 > lonsin[-1]:
        raise ValueError('lon0 outside of range of lonsin')
    i0 = np.argmin(np.fabs(lonsin-lon0))
    i0_shift = len(lonsin)-i0
    if ma.isMA(datain):
        dataout  = ma.zeros(datain.shape,datain.dtype)
    else:
        dataout  = np.zeros(datain.shape,datain.dtype)
    if ma.isMA(lonsin):
        lonsout = ma.zeros(lonsin.shape,lonsin.dtype)
    else:
        lonsout = np.zeros(lonsin.shape,lonsin.dtype)
    if start:
        lonsout[0:i0_shift] = lonsin[i0:]
    else:
        lonsout[0:i0_shift] = lonsin[i0:]-cyclic
    dataout[...,0:i0_shift] = datain[...,i0:]
    if start:
        lonsout[i0_shift:] = lonsin[start_idx:i0+start_idx]+cyclic
    else:
        lonsout[i0_shift:] = lonsin[start_idx:i0+start_idx]
    dataout[...,i0_shift:] = datain[...,start_idx:i0+start_idx]
    return dataout,lonsout

    def test_topo(self):
        file = self.f

        # basic case
        cstr_list = ('time|i0 zc|ZP|2500 yc|i5 xc|:', \
                'time|i0 zc|ZP|2500m yc|i5 xc|:', \
                'time|i0 zc|ZP|1500m yc|i5 xc|:', \
                'time|i0 zc|ZP|1000,1500m yc|i5.5 xc|:')
        results = ((21,), (21,), (21,), (2,21))

        for (cstr, res) in zip(cstr_list, results):
            if verbose: print cstr
            print "in test_topo"
            xsel = Nio.inp2xsel(file, 'PT', cstr)
            pt = file.variables['PT'][cstr]
            #pt = file.variables['ZP'][:]
            if verbose: print pt.shape
            if verbose: 
                if ma.isMA(pt):
                    print N.asarray(pt.filled())
                else:
                    print pt
            assert_equal(pt.shape, res)

        # ERROR:
        #cstr = 'xc|10k yc|i5.5:8:0.5i zc|ZP|2.5,3.5 time|i0:6:3'
        #if verbose: print cstr
        #pt = file.variables['PT'][cstr]
        #if verbose: print pt.shape

        file.close()

 def set_data(self, x, y, A):
     x = np.asarray(x,np.float32)
     y = np.asarray(y,np.float32)
     if not ma.isMA(A):
         A = np.asarray(A)
     if len(x.shape) != 1 or len(y.shape) != 1\
        or A.shape[0:2] != (y.shape[0], x.shape[0]):
         raise TypeError("Axes don't match array shape")
     if len(A.shape) not in [2, 3]:
         raise TypeError("Can only plot 2D or 3D data")
     if len(A.shape) == 3 and A.shape[2] not in [1, 3, 4]:
         raise TypeError("3D arrays must have three (RGB) or four (RGBA) color components")
     if len(A.shape) == 3 and A.shape[2] == 1:
         A.shape = A.shape[0:2]
     if len(A.shape) == 2:
         if A.dtype != np.uint8:
             A = (self.cmap(self.norm(A))*255).astype(np.uint8)
         else:
             A = np.repeat(A[:,:,np.newaxis], 4, 2)
             A[:,:,3] = 255
     else:
         if A.dtype != np.uint8:
             A = (255*A).astype(np.uint8)
         if A.shape[2] == 3:
             B = zeros(tuple(list(A.shape[0:2]) + [4]), np.uint8)
             B[:,:,0:3] = A
             B[:,:,3] = 255
             A = B
     self._A = A
     self._Ax = x
     self._Ay = y
     self._imcache = None

    def make_gene_map_2(self):
        """
        The method that takes the attributes from the array and uses
        them to create a gene map for the array.
        The gene map is a dictionary which has a binary string as a key.
        The binary string is created by creating a binary bit string of 
        an appropriate length.
        The length is calculated 
        """
        count = 0
        self.iterator_one = itertools.product(range(self.time_len), range(self.lat_len), range(self.lon_len))
        ### Assign a binary string a location and a value from the data
        print "\n"
        print "Creating gene-map dictionary... \n"
        print "Assigning valid locations to binary strings! \n"
        for x_valid in self.iterator_one:
            binary_string = bin(count)[2:]
            while len(binary_string) < self.string_length:  # removed minus one (-1) NB
                binary_string = "0" + binary_string
            self.gene_map[binary_string] = {}
            if ma.is_masked(self.array[x_valid]):
                pass
            else:
                self.gene_map[binary_string]["coordinate"] = tuple(x_valid)
                self.gene_map[binary_string]["value"] = self.array[x_valid]
                self.location_dict[x_valid[1:3]] = []
                self.location_dict_stdevs[x_valid[1:3]] = 0
                count += 1
        self.last_valid_binary_string = binary_string
        binary_string_old = binary_string
        not_valid_first = int(binary_string, 2) + 1
        not_valid_last = int("1" * (self.string_length), 2)  # added minus one just for nonmasked version NB
        self.count = count

        if self.count == self.count_non_masked:
            print "The counter corresponds with the non-masked count! \n"
        ### Pad the dictionary to give binary strings some value
        print "Assigning left over binary strings to non-existant locations! \n"
        self.iterator_two = itertools.product(range(self.time_len), range(self.lat_len), range(self.lon_len))
        # ~ for x_not_valid in range(not_valid_first, not_valid_last+1):
        count_2 = not_valid_first
        for x_not_valid in self.iterator_two:
            # ~ binary_string = bin(x_not_valid)[2:] # DOES IT NEED TO BE PADDED
            binary_string = bin(count_2)[2:]
            while len(binary_string) < self.string_length:  # removed minus one (-1) NB
                binary_string = "0" + binary_string
            self.gene_map[binary_string] = {}
            self.gene_map[binary_string]["coordinate"] = (999, 999, 999)  # x_not_valid
            self.gene_map[binary_string]["value"] = 1e09  # self.array[x_valid]
            if count_2 == not_valid_last:
                break
            else:
                count_2 += 1
        print "There are %d valid locations. \n" % count
        print "The last binary string is: ", binary_string
        print "The last binary string assigned to a valid locations is :", binary_string_old
        print "The length of binary string is: ", self.string_length
        print "The non-valid locations fall between %d and %d. \n" % (not_valid_first, not_valid_last)
        print "Is the array masked?: \n", ma.isMA(self.array)
        print "The gene-map has been created! \n"

def test_masked_unweighted():
    data_in = ma.ones((10,), dtype=[('x', float), ('y', float)])
    data_out = downsample(data_in, 2)
    assert ma.isMA(data_out)
    assert np.array_equal(data_out, data_in[:5])
    data_in['x'].mask[2] = True
    data_in.mask[7] = (True, True)
    data_out = downsample(data_in, 2)
    assert np.array_equal(data_out, data_in[:5])

def test_extrapolate():
    data = np.empty((10,), dtype=[('x', float), ('y', float)])
    data['x'] = np.arange(10.)
    data['y'] = np.arange(10.)
    x2 = np.arange(-1,11)
    result = resample(data, 'x', x2, 'y')
    assert ma.isMA(result)
    assert result['y'].mask[0]
    assert result['y'].mask[-1]
    assert np.array_equal(result['y'][1:-1], x2[1:-1])

def test_one_masked():
    data1 = ma.ones((10,), dtype=[('f', float), ('w', float)])
    data2 = np.ones((10,), dtype=[('f', float), ('w', float)])
    data1.mask = False
    data1['f'].mask[2] = True
    result = accumulate(data1_in=data1, data2_in=data2, add='f', weight='w')
    assert not ma.isMA(result), 'Result should not be masked.'
    assert np.all(result['f'] == 1), 'Incorrect addition result.'
    assert np.array_equal(result['w'][1:4], (2, 1, 2)),\
        'Mask not used correctly.'

def test_masked_weighted():
    data_in = ma.ones((10,), dtype=[('x', float), ('y', float)])
    data_out = downsample(data_in, 2, weight='y')
    assert ma.isMA(data_out)
    assert np.all(data_out['x'] == 1.)
    assert np.all(data_out['y'] == 2.)
    data_in['x'].mask[2] = True
    data_in.mask[7] = (True, True)
    data_out = downsample(data_in, 2, weight='y')
    assert np.all(data_out['x'] == 1.)
    assert np.all(data_out['y'] == (2., 1., 2., 1., 2.))

def addcyclic(arrin,lonsin):
    """
    ``arrout, lonsout = addcyclic(arrin, lonsin)``
    adds cyclic (wraparound) point in longitude to ``arrin`` and ``lonsin``.
    """
    nlats = arrin.shape[0]
    nlons = arrin.shape[1]
    if ma.isMA(arrin):
        arrout  = ma.zeros((nlats,nlons+1),arrin.dtype)
    else:
        arrout  = numpy.zeros((nlats,nlons+1),arrin.dtype)
    arrout[:,0:nlons] = arrin[:,:]
    arrout[:,nlons] = arrin[:,0]
    if ma.isMA(lonsin):
        lonsout = ma.zeros(nlons+1,lonsin.dtype)
    else:
        lonsout = numpy.zeros(nlons+1,lonsin.dtype)
    lonsout[0:nlons] = lonsin[:]
    lonsout[nlons]  = lonsin[-1] + lonsin[1]-lonsin[0]
    return arrout,lonsout

def test_propagated_array_mask():
    wlen = np.arange(10)
    flux = ma.ones((10,))
    flux.mask = False
    flux[2] = ma.masked
    result = redshift(z_in=0, z_out=1, rules=[
        {'name': 'wlen', 'exponent': +1, 'array_in': wlen},
        {'name': 'flux', 'exponent': -1, 'array_in': flux}])
    assert ma.isMA(result)
    assert not result['wlen'].mask[2], 'Input mask not propagated.'
    assert result['flux'].mask[2], 'Input mask not propagated.'
    assert not result['flux'].mask[3], 'Input mask not propagated.'

def test_propagated_data_mask():
    data_in = ma.ones((10,), dtype=[
        ('wlen', float), ('flux', float), ('extra', int)])
    data_in['wlen'][1] = ma.masked
    data_in['extra'][2] = ma.masked
    result = redshift(z_in=0, z_out=1, data_in=data_in, rules=[
        {'name': 'wlen', 'exponent': +1},
        {'name': 'flux', 'exponent': -1}])
    assert ma.isMA(result)
    assert not result['wlen'].mask[0], 'Input mask not propagated.'
    assert not result['flux'].mask[0], 'Input mask not propagated.'
    assert result['wlen'].mask[1], 'Input mask not propagated.'
    assert result['extra'].mask[2], 'Input mask not propagated.'

    def __init__(self, var, name=None):
        """Create a VarInfo object.

        Arguments:
        var: numpy or numpy.ma array"""

        # Compute all necessary statistics in initialization, so that we don't
        # have to hold onto the variable in memory for later use (in case the
        # variable consumes a lot of memory).
        if not ma.isMA(var):
            var = ma.array(var)
        self._compute_stats(var)
        self.name = name

def test_both_masked():
    data1 = ma.ones((10,), dtype=[('f', float), ('w', float), ('i', int)])
    data2 = ma.ones((10,), dtype=[('f', float), ('w', float), ('i', int)])
    data1.mask = False
    data1['f'].mask[2:4] = True
    data2.mask = False
    data2['f'].mask[3:5] = True
    result = accumulate(data1_in=data1, data2_in=data2,
                        add='f', weight='w', join='i')
    assert not ma.isMA(result), 'Result should not be masked.'
    valid = result['w'] != 0
    assert np.all(result['f'][valid] == 1), 'Incorrect addition result.'
    assert np.array_equal(result['w'][1:6], (2, 1, 0, 1, 2)),\
        'Mask not used correctly.'

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 fill_gaps_in_2d_transect_once(Z):
    was_masked = ma.isMA(Z)

    # Ensure array has NaNs, not mask
    Z = ma.filled(Z, np.nan).copy()

    fill_values = np.nanmean(
        np.dstack((np.roll(Z, 1, axis=1), np.roll(Z, -1, axis=1))), axis=2)

    inds_to_fill = np.logical_and(np.isnan(Z), ~np.isnan(fill_values))

    Z[inds_to_fill] = fill_values[inds_to_fill]

    if was_masked:
        Z = ma.masked_invalid(Z)

    return Z

def _quantize(data,least_significant_digit):
    """
quantize data to improve compression. data is quantized using
around(scale*data)/scale, where scale is 2**bits, and bits is determined
from the least_significant_digit. For example, if
least_significant_digit=1, bits will be 4.
    """
    precision = pow(10.,-least_significant_digit)
    exp = np.log10(precision)
    if exp < 0:
        exp = int(np.floor(exp))
    else:
        exp = int(np.ceil(exp))
    bits = np.ceil(np.log2(pow(10.,-exp)))
    scale = pow(2.,bits)
    datout = np.around(scale*data)/scale
    if ma.isMA(datout):
        datout.set_fill_value(data.fill_value)
        return datout
    else:
        return datout

def bin_2d_transect(x, y, Z, x_out, y_out):
    """Bin transect Z(x, y), where x can be irregular

    Inputs
    ------
    x, y : 1D arrays
        x can be irregular, y cannot
    Z : 2D array
        Data at each point x, y. May be masked array
    x_out, y_out : 1D arrays
        Edges of grid on which to bin Z

    Returns
    -------
    Z_out : 2D array
        Shape (len(x_out) - 1, len(y_out) - 1)
    """
    if Z.ndim == 1:
        Z = Z[np.newaxis, :]

    # Preallocate result
    Nx, Ny = x_out.size - 1, y_out.size - 1
    Z_out = np.full((Nx, Ny), np.nan)

    filterwarnings('ignore', '.*Mean of empty slice*.')

    # Using loop for simplicity
    for i, j in np.ndindex(Nx, Ny):
        in_x_bin = np.logical_and(x > x_out[i], x < x_out[i + 1])
        in_y_bin = np.logical_and(y > y_out[j], y < y_out[j + 1])

        Z_in_bin = Z[in_y_bin, in_x_bin]
        Z_out[i, j] = np.nanmean(ma.filled(Z_in_bin, np.nan))

    if ma.isMA(Z):
        Z_out = ma.masked_invalid(Z_out)

    return Z_out

    def set_data(self, A, shape=None):
        """
        Set the image array

        ACCEPTS: numpy/PIL Image A"""
        # check if data is PIL Image without importing Image
        if hasattr(A,'getpixel'):
            self._A = pil_to_array(A)
        elif ma.isMA(A):
            self._A = A
        else:
            self._A = np.asarray(A) # assume array

        if self._A.dtype != np.uint8 and not np.can_cast(self._A.dtype, np.float):
            raise TypeError("Image data can not convert to float")

        if (self._A.ndim not in (2, 3) or
            (self._A.ndim == 3 and self._A.shape[-1] not in (3, 4))):
            raise TypeError("Invalid dimensions for image data")

        self._imcache =None
        self._rgbacache = None
        self._oldxslice = None
        self._oldyslice = None

def shiftgrid(lon0,datain,lonsin,start=True,cyclic=360.0):

    import numpy.ma as ma
    """
    Shift global lat/lon grid east or west.
    copied directly from mpl_toolkits v1.0.2 by mjh

    .. tabularcolumns:: |l|L|

    ==============   ====================================================
    Arguments        Description
    ==============   ====================================================
    lon0             starting longitude for shifted grid
                     (ending longitude if start=False). lon0 must be on
                     input grid (within the range of lonsin).
    datain           original data.
    lonsin           original longitudes.
    ==============   ====================================================

    .. tabularcolumns:: |l|L|

    ==============   ====================================================
    Keywords         Description
    ==============   ====================================================
    start            if True, lon0 represents the starting longitude
                     of the new grid. if False, lon0 is the ending
                     longitude. Default True.
    cyclic           width of periodic domain (default 360)
    ==============   ====================================================

    returns ``dataout,lonsout`` (data and longitudes on shifted grid).
    """
    if numpy.fabs(lonsin[-1]-lonsin[0]-cyclic) > 1.e-4:
        # Use all data instead of raise ValueError, 'cyclic point not included'
        start_idx = 0
    else:
        # If cyclic, remove the duplicate point
        start_idx = 1
    if lon0 < lonsin[0] or lon0 > lonsin[-1]:
        msg = 'lon0 outside of range of lonsin %(l0)4.1f %(st)4.1f %(ed)4.1f'%{'l0':lon0,'st':lonsin[0],'ed':lonsin[-1]}
        raise ValueError(msg)
    i0 = numpy.argmin(numpy.fabs(lonsin-lon0))
    i0_shift = len(lonsin)-i0
    if ma.isMA(datain):
        dataout  = ma.zeros(datain.shape,datain.dtype)
    else:
        dataout  = numpy.zeros(datain.shape,datain.dtype)
    if ma.isMA(lonsin):
        lonsout = ma.zeros(lonsin.shape,lonsin.dtype)
    else:
        lonsout = numpy.zeros(lonsin.shape,lonsin.dtype)
    if start:
        lonsout[0:i0_shift] = lonsin[i0:]
    else:
        lonsout[0:i0_shift] = lonsin[i0:]-cyclic
    dataout[:,0:i0_shift] = datain[:,i0:]
    if start:
        lonsout[i0_shift:] = lonsin[start_idx:i0+start_idx]+cyclic
    else:
        lonsout[i0_shift:] = lonsin[start_idx:i0+start_idx]
    dataout[:,i0_shift:] = datain[:,start_idx:i0+start_idx]
    return dataout,lonsout