def spike_flag(data,masked,freq,percent):
    """
    Flags out RFI spikes using a 11 bin filter
    Can be used with either time or freq
    percent is a percentage level cut (100 would be twice the 11 bin average)
    Needs to be applied to masked data.
    """
    new_mask = np.zeros(len(data))
    new_array = ma.array(data,mask=masked)
    new_comp = ma.compressed(new_array)
    freq_array = ma.array(freq,mask=masked)
    new_freq = ma.compressed(freq_array)
    for i in range(0,len(data)):
        if masked[i]==1.0:
            new_mask[i] = 1.0

    for i in range(5,len(new_comp)-5):
        group = new_comp[i-5]+new_comp[i-4]+new_comp[i-3]+new_comp[i-2]+new_comp[i-1]+new_comp[i]+new_comp[i+1]+new_comp[i+2]+new_comp[i+3]+new_comp[i+4]+new_comp[i+5]
        mean_group = group/11.
        if new_comp[i]/mean_group>=(1+percent/100.):
            comp_freq = new_freq[i]
            for j in range(0,len(freq)):
                if freq[j]==comp_freq:
                    index=j
            new_mask[index]= 1.0
        elif new_comp[i]/mean_group<=1/(1+percent/100.):
            comp_freq = new_freq[i]
            for j in range(0,len(freq)):
                if freq[j]==comp_freq:
                    index=j
            new_mask[index]= 1.0
   
    return new_mask

    def set_gen3data(self, Apath, Ppath, Xdet, Zdet, 
            omega, mask_cutoff, toprowsmask=0, polygonmaskfile=None):

        self.omega = omega
        self.wbyv = omega/self.v

        # read the data
        self.Amp = bp.io.read_frame(Apath)
        self.Pha = bp.io.read_frame(Ppath)
        self.amshape = self.Amp.shape

        # set the mask
        self.mask = np.ones_like(self.Amp, dtype='bool')
        self.mask[self.Amp > mask_cutoff*np.max(self.Amp)] = False
        # mask top rows
        if toprowsmask > 0:
            toprowsmask = toprowsmask -1 
            self.mask[0:toprowsmask,:] = True

        # mask all data arrays
        self.Amp = ma.compressed((ma.masked_array(self.Amp, self.mask)))
        self.Pha = ma.compressed(ma.masked_array(self.Pha, self.mask))

        # masked detector positions
        self.dIdx = np.arange(len(Xdet))
        self.dIdx = ma.compressed(ma.masked_array(self.dIdx.reshape(self.amshape), self.mask))
        self.Xdet = ma.compressed(ma.masked_array(Xdet.reshape(self.amshape), self.mask))
        self.Zdet = ma.compressed(ma.masked_array(Zdet.reshape(self.amshape), self.mask))

    def set_gen3data(self, Apaths, Ppaths, Xdet, Zdet, omega, mask=None):

        self.omega = omega
        self.wbyv = omega/self.v
        self.amshape = mask.shape
        self.mask = mask

        for Apath, Ppath in zip(Apaths, Ppaths):
            # read the data
            Amp = bp.io.read_frame(Apath)
            Pha = bp.io.read_frame(Ppath)

            # mask all data arrays
            Amp = ma.compressed(ma.masked_array(Amp, self.mask))
            Pha = ma.compressed(ma.masked_array(Pha, self.mask))

            try:
                Amps = np.vstack((Amps, Amp))
                Phas = np.vstack((Phas, Pha))
            except:
                Amps = Amp
                Phas = Pha
        self.Amps = Amps
        self.Phas = Phas

        # masked detector positions
        self.dIdx = np.arange(len(Xdet))
        self.dIdx = ma.compressed(ma.masked_array(self.dIdx.reshape(self.amshape), self.mask))
        self.Xdet = ma.compressed(ma.masked_array(Xdet.reshape(self.amshape), self.mask))
        self.Zdet = ma.compressed(ma.masked_array(Zdet.reshape(self.amshape), self.mask))

    def calculateMeans(self):
        self.synHist = ma.masked_values(self.synHist, -9999.)
        self.synHistMean = ma.mean(self.synHist, axis=0)
        self.medSynHist = ma.median(self.synHist, axis=0)

        self.synHistUpper = percentile(ma.compressed(self.synHist),
                                       per=95, axis=0)
        self.synHistLower = percentile(ma.compressed(self.synHist),
                                       per=5, axis=0)

 def Fill2ThetaMap(data,TA,image):
     import numpy.ma as ma
     Zmin = data['Zmin']
     Zmax = data['Zmax']
     tax,tay = TA    # 2-theta & yaxis
     taz = ma.masked_outside(image.flatten()-Zmin,0,Zmax-Zmin)
     tam = ma.getmask(taz)
     tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
     tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
     taz = ma.compressed(ma.array(taz.flatten(),mask=tam))
     del(tam)
     return tax,tay,taz

def azimuthalAverage(image, center=None, maskval=0):
    """
    calculate the azimuthally averaged radial profile.

    image - 2D image
    center - [x,y] pixel coordinates used as the center. the default is
             None which then uses the center of the image (including
             fractional pixels).
    maskval - threshold value for including data in the profile
    """

    # calculate the indices from the image
    y, x = np.indices(image.shape)

    # default to image center if no center given
    if not center:
        center = np.array([(x.max() - x.min()) / 2.0,
                           (x.max() - x.min()) / 2.0])

    r = np.hypot(x - center[0], y - center[1])

    # get sorted radii and sort image accordingly
    ind = np.argsort(r.flat)
    i_sorted = image.flat[ind]

    # for FP data we need to at least mask out data at
    # 0 or less so the gaps get ignored.
    # also want to mask out area outside of aperture
    # so use given maskval to do that.
    i_ma = ma.masked_less_equal(i_sorted, maskval)
    mask = ma.getmask(i_ma)

    # remove masked data points from further analysis
    r_sorted = ma.compressed(ma.array(r.flat[ind], mask=mask))
    i_mask = ma.compressed(i_ma)

    # get the integer part of the radii (bin size = 1)
    r_int = r_sorted.astype(int)

    # find all pixels that fall within each radial bin.
    deltar = r_int[1:] - r_int[:-1]  # assumes all radii represented
    rind = np.where(deltar)[0]       # location of changed radius
    nr_tot = rind[1:] - rind[:-1]    # total number of points in radius bin

    # cumulative sum to figure out sums for each radius bin
    csim = ma.cumsum(i_mask, dtype=float)
    tbin = csim[rind[1:]] - csim[rind[:-1]]

    # calculate and return profile of mean within each bin
    radial_prof = tbin / nr_tot

    return radial_prof

def loessmin_diff((file1, file2, filex1, filex2, cutoff, toprowsmask, span)):
    from bopy.io import read_frame
    from bopy.utils import mask2_maxcutoff 
    from bopy.utils.rfuncs import loess2d
    d1 = read_frame(file1)
    d2 = read_frame(file2)
    dmask = mask2_maxcutoff(read_frame(filex1),
            read_frame(filex2),
            cutoff, toprowsmask=toprowsmask)
    d1 = loess2d(d1, span=span)
    d1 = np.min(ma.compressed(ma.masked_array(d1, mask=dmask)))
    d2 = loess2d(d2, span=span)
    d2 = np.min(ma.compressed(ma.masked_array(d2, mask=dmask)))
    return d1 - d2

def band2txt(band, outfile):
    """Accepts numpy rasterand writes to specified text file on disk."""
    if ma.isMaskedArray(band) is True:
        outraster = ma.compressed(band)
    else:
        outraster = band
    np.savetxt(outfile, outraster, fmt="%f")

def EdgeFinder(image,data):
    '''this makes list of all x,y where I>edgeMin suitable for an ellipse search?
    Not currently used but might be useful in future?
    '''
    import numpy.ma as ma
    Nx,Ny = data['size']
    pixelSize = data['pixelSize']
    edgemin = data['edgemin']
    scalex = pixelSize[0]/1000.
    scaley = pixelSize[1]/1000.    
    tay,tax = np.mgrid[0:Nx,0:Ny]
    tax = np.asfarray(tax*scalex,dtype=np.float32)
    tay = np.asfarray(tay*scaley,dtype=np.float32)
    tam = ma.getmask(ma.masked_less(image.flatten(),edgemin))
    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
    return zip(tax,tay)

def gain_calc(data,masked,gsm,K0):
    """
    Calculates the gain using least squares for a single frequency.
    Inputs:
    data - single freq array of data
    gsm  - similar single freq array of gsm data
    masked - mask for data
    K0 - preliminary guess for gain
    """
    fK = lambda K,d,g: K*d-g

    d0_array = ma.array(data,mask=masked)
    d0 = ma.compressed(d0_array)
    g0_array = ma.array(gsm,mask=masked)
    g0 = ma.compressed(g0_array)
    Kf = opt.leastsq(fK,K0,args=(d0,g0),maxfev=100000)
    
    return Kf[0]

def peak_finder(arr, smoothing=5, ddy_thresh=-300, dy0_thresh=5):
    array        = medfilt(arr, smoothing) 
    x            = arange(len(array))
    kernel       = [4, 0, -4]
    dY           = convolve(array, kernel, 'same') 
    ddy          = convolve(dY, kernel, 'same')
    falloff      = -15000*exp(-0.003*x) #This has to be worked on
    masked_array = ma.masked_where(logical_or(ddy>falloff+ddy_thresh, abs(dY) > dy0_thresh) , arr) 
    x_masked     = ma.array(x, mask=masked_array.mask)
    return ma.compressed(x_masked)

    def predict(self, X):
        # Handle missing values
        X = ma.masked_invalid(X)
        mask = X.mask
        dX = ma.compressed(X).reshape(-1, 1)
        dZ = self.model.predict(dX)
        Z = np.array([np.nan for i in xrange(X.shape[0])])
        Z[~mask] = dZ
        Z = ma.masked_invalid(Z)

        return Z * self.step

    def calculateMeans(self, synMean, synMin, synMed, synMax, synMinCP):
        synMean = ma.masked_values(synMean, -9999.)
        synMin = ma.masked_values(synMin, -9999.)
        synMed = ma.masked_values(synMed, -9999.)
        synMax = ma.masked_values(synMax, -9999.)

        self.synMean = ma.mean(synMean, axis=0)
        self.synMed = ma.mean(synMed, axis=0)
        self.synMin = ma.mean(synMin, axis=0)
        self.synMax = ma.mean(synMax, axis=0)

        self.synMeanUpper = percentile(ma.compressed(synMean), per=95, axis=0)
        self.synMeanLower = percentile(ma.compressed(synMean), per=5, axis=0)
        self.synMinUpper = percentile(ma.compressed(synMin), per=95, axis=0)
        self.synMinLower = percentile(ma.compressed(synMin), per=5, axis=0)

        self.synMinCPDist = np.mean(synMinCP, axis=0)
        self.synMinCPLower = percentile(synMinCP, per=5, axis=0)
        self.synMinCPUpper = percentile(synMinCP, per=95, axis=0)
        
        r = list(np.random.uniform(high=synMean.shape[0], size=3).astype(int))
        self.synRandomMinima = synMean[r, :, :]

def Fill2ThetaAzimuthMap(masks,TA,tam,image):
    'Needs a doc string'
    Zlim = masks['Thresholds'][1]
    rings = masks['Rings']
    arcs = masks['Arcs']
    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]),ma.getdata(TA[2])))    #azimuth, 2-theta, dist
    tax,tay,tad = np.dsplit(TA,3)    #azimuth, 2-theta, dist**2/d0**2
    for tth,thick in rings:
        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
    for tth,azm,thick in arcs:
        tamt = ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))
        tama = ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1]))
        tam = ma.mask_or(tam.flatten(),tamt*tama)
    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
    tabs = np.ones_like(taz)
    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))   #azimuth
    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))   #2-theta
    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))   #intensity
    tad = ma.compressed(ma.array(tad.flatten(),mask=tam))   #dist**2/d0**2
    tabs = ma.compressed(ma.array(tabs.flatten(),mask=tam)) #ones - later used for absorption corr.
    return tax,tay,taz,tad,tabs

def rat_fore(data,masked,freq,minf,maxf,n,m):
    """
    Calculates a rational function fit for the data.
    Inputs:
    data - single frequency dependent spectrum
    masked - corresponding mask
    freq - corresponding frequency array
    minf, maxf - min and max freq if want to truncate range
    n - index of numerator polynomial fitting
    m - index of denominator polynomial fitting
    (4>= n,m >=1)
    Output:
    dfit - polynomial fit spectrum
    fit_params - parameters for the fit
    """ 
    data_array = ma.array(data,mask=masked)
    data_comp = ma.compressed(data_array)
    freq_array = ma.array(freq,mask=masked)
    freq_comp = ma.compressed(freq_array)
 
    min_ind = 0
    max_ind = -1
    if minf>freq_comp[0]:
        min_ind = np.where(freq_comp<=minf)[0][-1]
    if maxf<freq_comp[-1]:
        max_ind = np.where(freq_comp<=maxf)[0][-1]
    mid_ind = min_ind+(max_ind-min_ind)/2 
 
    log_data = np.log10(data_comp[min_ind:max_ind])
    log_freq = np.log10(freq_comp[min_ind:max_ind]/freq_comp[mid_ind])
    
    p0 = np.ones(n+m+2)
    fitfunc = rational_fit(n,m)
    errfunc = lambda p,x,y: fitfunc(p,x)-y
    pf,success = opt.leastsq(errfunc,p0[:],args=(log_freq,log_data))
    dfit = 10**(fitfunc(pf,np.log10(freq/freq_comp[mid_ind])))

    return dfit,pf

def rebin(data,masked,freq,binscale):
    """
    Rebins data to coarser frequency resolution.
    Assumes that the input is the data after flagging, mask,
    corresponding freq array and a binscale (number of bins to merge).
    Output is rebinned data with corresponding freq, mask arrays.
    """

    if binscale > 1:
        new_data = np.zeros(len(data)/binscale)
        new_mask = np.zeros(len(data)/binscale)
        new_freq = np.zeros(len(data)/binscale)
        f=0
        for f in range(0, len(new_data)-1):
            if len(masked[f*binscale:(f+1)*binscale])==sum(masked[f*binscale:(f+1)*binscale]):
                new_data[f] = 1.0
            else: 
                test_data = ma.array(data[f*binscale:(f+1)*binscale],mask=masked[f*binscale:(f+1)*binscale])
                test_data_con = ma.compressed(test_data)
                new_data[f] = ma.mean(test_data_con)
            if sum(masked[f*binscale:(f+1)*binscale])>=binscale/2.:
                new_mask[f] = 1.0
            new_freq[f] = ma.mean(freq[f*binscale:(f+1)*binscale])
        if len(masked[(f+1)*binscale:-1])==sum(masked[(f+1)*binscale:-1]):
            new_data[-1] = 1.0
        else:
            test_data = ma.array(data[(f+1)*binscale:-1],mask=masked[(f+1)*binscale:-1])
            test_data_con = ma.compressed(test_data) 
            new_data[-1] = ma.mean(test_data_con)
        if sum(masked[(f+1)*binscale:-1])>=1.:
            new_mask[-1] = 1.0
        new_freq[-1] = ma.mean(freq[(f+1)*binscale:-1])
    elif binscale == 1:
        new_data = data
        new_mask = masked
        new_freq = freq
    
    return new_data,new_mask,new_freq

    def set_gen3data(self, Apath, Ppath, Xdet, Zdet, 
            omega, mask=None):

        self.omega = omega
        self.wbyv = omega/self.v

        # read the data
        self.Amp = bp.io.read_frame(Apath)
        self.Pha = bp.io.read_frame(Ppath)
        self.amshape = self.Amp.shape

        # set the mask
        self.mask = mask

        # mask all data arrays
        self.Amp = ma.compressed((ma.masked_array(self.Amp, self.mask)))
        self.Pha = ma.compressed(ma.masked_array(self.Pha, self.mask))

        # masked detector positions
        self.dIdx = np.arange(len(Xdet))
        self.dIdx = ma.compressed(ma.masked_array(self.dIdx.reshape(self.amshape), self.mask))
        self.Xdet = ma.compressed(ma.masked_array(Xdet.reshape(self.amshape), self.mask))
        self.Zdet = ma.compressed(ma.masked_array(Zdet.reshape(self.amshape), self.mask))

    def _get_bounds(item_list, wrapped_coords=False):
        """
        Return a tuple containing the first and secound bound in the list.
            * For 0 values it returns (None, None).
            * For 1 value it returns that value twice
            * For 2 values it returns those (low, high) for unwrapped
              co-ordinates or in the order given for wrapped
            * For Multiple values unwrapped, this returns min and max
            * For Multiple values wrapped co-ordinates, this
              returns the value around the largest gap in values

        :param list item_list: List of comparable data items
        :param wrapped_coords: is this a coordinate which wraps at 360 back to
                               0, e.g. longitude
        :return: Tuple of (first bound, second bound for values in the list)
        """
        items = ma.compressed(item_list)
        if len(items) < 1:
            return None, None

        if len(items) == 1:
            return items[0], items[0]

        if wrapped_coords:
            if len(items) is 2:
                first_bound_index = 0
                second_bound_index = 1
            else:
                # find the largest angle between closest points and exclude
                # this from the bounding box ensuring that this includes going
                # across the zero line
                items = sorted(items)
                first_bound_index = 0
                second_bound_index = len(items) - 1
                max_diff = (items[first_bound_index] -
                            items[second_bound_index]) % 360
                for i in range(1, len(items)):
                    diff = (items[i] - items[i-1]) % 360
                    if diff > max_diff:
                        max_diff = diff
                        first_bound_index = i
                        second_bound_index = i-1

            first_bound = items[first_bound_index]
            second_bound = items[second_bound_index]
        else:
            first_bound = min(items)
            second_bound = max(items)

        return float(first_bound), float(second_bound)

def poly_fore(data,masked,freq,minf,maxf,n,std):
    """
    Calculates a polynomial fit for data.
    Inputs:
    data - single frequency dependent spectrum
    masked - corresponding mask
    freq - corresponding frequency array
    minf, maxf - min and max freq if want to truncate range 
    n - index of polynomial fitting
    std - 1/weights, set to ones if not needed.
    Output:
    dfit - polynomial fit spectrum
    fit_params - parameters for the fit
    """
    
    data_array = ma.array(data,mask=masked)
    data_comp = ma.compressed(data_array)
    freq_array = ma.array(freq,mask=masked)
    freq_comp = ma.compressed(freq_array)
    std_comp = ma.compressed(ma.array(std, mask=masked))

    min_ind = 0
    max_ind = -1
    if minf>freq_comp[0]:
        min_ind = np.where(freq_comp<=minf)[0][-1]
    if maxf<freq_comp[-1]:
        max_ind = np.where(freq_comp<=maxf)[0][-1]
    mid_ind = min_ind+(max_ind-min_ind)/2

    log_data = np.log10(data_comp[min_ind:max_ind])
    log_freq = np.log10(freq_comp[min_ind:max_ind]/freq_comp[mid_ind])
    weights = 1/std_comp[min_ind:max_ind]
    
    fit_params = poly.polyfit(log_freq,log_data,n,w=weights)
    dfit = 10**(poly.polyval(np.log10(freq/freq_comp[mid_ind]),fit_params))

    return dfit, fit_params

def rasterHistogram(raster_matrix):
    '''Accepts matrix and generates histogram'''
     
    # Line above is function docstring
    # Flatten 2d-matrix
    flat_raster = ma.compressed(raster_matrix)
     
    # Setup the plot (see matplotlib.sourceforge.net)
    fig = plt.figure(figsize=(8,11))
    ax = fig.add_subplot(1,1,1)
     
    # Plot histogram
    ax.hist(flat_raster, 10, normed=0, histtype='bar',
            align='mid')
    # Show the plot on screen
    plt.show()