def __init__(self, cube, wcs=None, mask=None, sigma=None, empty_channel=0,
                 keep_threshold_mask=True, distance=None, galaxy_props={}):
        super(BubbleFinder, self).__init__()

        if not isinstance(cube, SpectralCube):
            if wcs is None:
                raise TypeError("When cube is not a SpectralCube, wcs must be"
                                " given.")
            cube = SpectralCube(cube, wcs)
            if mask is not None:
                cube = cube.with_mask(mask)

        self.cube = cube

        self.empty_channel = empty_channel

        if sigma is None:
            self.estimate_sigma()
        else:
            self.sigma = sigma

        self.keep_threshold_mask = keep_threshold_mask
        self._mask = None
        self.distance = distance
        self.galaxy_props = galaxy_props

def line_flux2(catalog, line_name='13co10',
               asgn=datadir + 'COHRS_all_asgn.fits',
               cubefile=datadir + 'GRS_13CO_all.fits'):

    flux = Column(np.zeros(len(catalog)),name=line_name)

    asgn = SpectralCube.read(asgn)
    linefile = SpectralCube.read(cubefile)

    previous_file=''
    fill_data=None
    previous_cube_file=''
    for idx, obj in enumerate(catalog):
        if obj['orig_file'] != previous_cube_file:
            print "Pulling line subcube for {0}".format(obj['orig_file'])
            subx1 = obj['orig_file'].split('_')[2]
            subx2 = obj['orig_file'].split('_')[3]
            subcube = linefile[:, :, int(subx1):int(subx2)]
            fill_cube_data = (subcube.filled_data[:].value)
            previous_cube_file = obj['orig_file']
        
        outtuple = sparse_mask(obj, asgn,
                               previous_file=previous_file,
                               fill_data=fill_data)
        previous_file, fill_data, zcld, ycld, xcld = outtuple
        if len(xcld)>0:
            flux[idx] = np.nansum(fill_cube_data[zcld, ycld, xcld])
    catalog.add_column(flux)
    return catalog

def load_and_reduce(filename, add_noise=False, rms_noise=0.001,
                    nsig=3):
    '''
    Load the cube in and derive the property arrays.
    '''

    if add_noise:
        if rms_noise is None:
            raise TypeError("Must specify value of rms noise.")

        cube, hdr = getdata(filename, header=True)

        from scipy.stats import norm
        cube += norm.rvs(0.0, rms_noise, cube.shape)

        sc = SpectralCube(data=cube, wcs=WCS(hdr))

        mask = LazyMask(np.isfinite, sc)
        sc = sc.with_mask(mask)

    else:
        sc = filename

    reduc = Mask_and_Moments(sc, scale=rms_noise)
    reduc.make_mask(mask=reduc.cube > nsig * reduc.scale)
    reduc.make_moments()
    reduc.make_moment_errors()

    return reduc.to_dict()

    def __init__(self, highres, lowres):
        super(MultiResObs, self).__init__()
        self.highres = SpectralCube.read(highres)
        self.lowres = SpectralCube.read(lowres)

        self.highres_convolved = None
        self.lowres_convolved = None

        self.lowbeam = self.lowres.beam
        self.highbeam = self.highres.beam

        self.combined_beam = self.lowbeam.convolve(self.highbeam)

def cubegen(ymin,ymax,xmin,xmax, deltaX=30):
	"""Generates a subcube of the specified dimensions from the .fits files,
	   for 12CO and 13CO. Returns the subcubes for 12CO and 13CO, respectively.

	   Argument format: "(ymin,ymax, xmin,xmax.)"
	   ^ These are the parameters of the desired subcubes."""

	cube12 = SpectralCube.read("paws-30m-12co10-23as-cube.fits")
	cube13 = SpectralCube.read("paws-30m-13co10-23as-cube.fits")
	subcube12 = cube12[:,ymin:ymax,xmin:xmax]
	subcube13 = cube13[:,ymin:ymax,xmin:xmax]

	return subcube12,subcube13

def reduce_and_save(filename, add_noise=False, rms_noise=0.001,
                    output_path="", cube_output=None,
                    nsig=3, slicewise_noise=True):
    '''
    Load the cube in and derive the property arrays.
    '''

    if add_noise:
        if rms_noise is None:
            raise TypeError("Must specify value of rms noise.")

        cube, hdr = getdata(filename, header=True)

        # Optionally scale noise by 1/10th of the 98th percentile in the cube
        if rms_noise == 'scaled':
            rms_noise = 0.1*np.percentile(cube[np.isfinite(cube)], 98)

        from scipy.stats import norm
        if not slicewise_noise:
            cube += norm.rvs(0.0, rms_noise, cube.shape)
        else:
            spec_shape = cube.shape[0]
            slice_shape = cube.shape[1:]
            for i in range(spec_shape):
                cube[i, :, :] += norm.rvs(0.0, rms_noise, slice_shape)

        sc = SpectralCube(data=cube, wcs=WCS(hdr))

        mask = LazyMask(np.isfinite, sc)
        sc = sc.with_mask(mask)

    else:
        sc = filename

    reduc = Mask_and_Moments(sc, scale=rms_noise)
    reduc.make_mask(mask=reduc.cube > nsig * reduc.scale)

    reduc.make_moments()
    reduc.make_moment_errors()

    # Remove .fits from filename
    save_name = filename.split("/")[-1][:-4]

    reduc.to_fits(output_path+save_name)

    # Save the noisy cube too
    if add_noise:
        if cube_output is None:
            reduc.cube.hdu.writeto(output_path+save_name)
        else:
            reduc.cube.hdu.writeto(cube_output+save_name)

def warp_ellipse_to_circle(cube, a, b, pa, stop_if_huge=True):
    '''
    Warp a SpectralCube such that the given ellipse is a circle int the
    warped frame.

    Since you should **NOT** be doing this with a large cube, we're going
    to assume that the given cube is a subcube centered in the middle of the
    cube.

    This requires a rotation, then scaling. The equivalent matrix is:
    [b cos PA    b sin PA]
    [-a sin PA   a cos PA ].

    '''

    if cube._is_huge:
        if stop_if_huge:
            raise Warning("The cube has the huge flag enabled. Disable "
                          "'stop_if_huge' if you would like to continue "
                          "anyways with the warp.")
        else:
            warn("The cube has the huge flag enabled. This may use a lot "
                 "of memory!")

    # Let NaNs be 0
    data = cube.with_fill_value(0.0).filled_data[:].value

    warped_array = []

    for i in range(cube.shape[0]):
        warped_array.append(nd.zoom(nd.rotate(data[i], np.rad2deg(-pa)),
                                    (1, a / b)))

    warped_array = np.array(warped_array)

    # We want to mask outside of the original bounds
    mask = np.ones(data.shape[1:])
    warp_mask = \
        np.isclose(nd.zoom(nd.rotate(mask, np.rad2deg(-pa)),
                           (1, a / b)), 1)

    # There's probably a clever way to transform the WCS, but all the
    # solutions appear to need pyast/starlink. The output of the wrap should
    # give a radius of b and the spectral dimension is unaffected.
    # Also this is hidden and users won't be able to use this weird cube
    # directly
    warped_cube = SpectralCube(warped_array * cube.unit, cube.wcs)
    warped_cube = warped_cube.with_mask(warp_mask)

    return warped_cube

def select_cloud(idxarray, cloudcat):
    for idx in idxarray:
        entry = cloudcat[idx]
        asgn = SpectralCube.read(cohrsdir+'FINALASGNS/'+
                                 entry['orig_file']+
                                 '_fasgn.fits')

        data = SpectralCube.read(cohrsdir+'DATA/'+
                                 entry['orig_file']+
                                 '.fits')
        mask = (asgn == entry['_idx'] *
                u.dimensionless_unscaled)
        cube = data.with_mask(mask)
        cube = cube.minimal_subcube()
        cube.write('cohrscld_{0}'.format(entry['_idx'])+'.fits',
                   overwrite=True)

def test_qglue():
    from spectral_cube import SpectralCube
    cube = SpectralCube.read(os.path.join(DATA, 'cube_3d.fits'))
    data = parse_data(cube, 'x')[0]
    assert data.label == 'x'
    data['STOKES I']
    assert data.shape == (2, 3, 4)

def FirstLook_Cepheus():
    print("Now NH3(1,1)")
    a_rms = [  0, 135, 290, 405, 505, 665]
    b_rms = [ 70, 245, 350, 455, 625, 740]
    index_rms=first_look.create_index( a_rms, b_rms)
    index_peak=np.arange(350,410)
    file_in='Cepheus/Cepheus_NH3_11.fits'
    # 1st order polynomial
    file_out=file_in.replace('.fits','_base1.fits')
    file_new=first_look.baseline( file_in, file_out, index_clean=index_rms, polyorder=1)
    first_look.peak_rms( file_new, index_rms=index_rms, index_peak=index_peak)
    print("Now NH3(2,2)")
    linelist = ['NH3_22','NH3_33','C2S','HC5N','HC7N_21_20','HC7N_22_21']
    vsys = -3.8*u.km/u.s
    throw = 2.0*u.km/u.s
    for line in linelist:
        file_in = 'Cepheus/Cepheus_{0}.fits'.format(line)
        s = SpectralCube.read(file_in)
        s = s.with_spectral_unit(u.km/u.s,velocity_convention='radio')
        a_rms = [s.closest_spectral_channel(vsys+3*throw),
                 s.closest_spectral_channel(vsys-throw)]
        b_rms = [s.closest_spectral_channel(vsys+throw),
                 s.closest_spectral_channel(vsys-3*throw)]
        index_peak = np.arange(s.closest_spectral_channel(vsys+3*u.km/u.s),
                              s.closest_spectral_channel(vsys-3*u.km/u.s))
        index_rms=first_look.create_index( a_rms, b_rms)

        file_out=file_in.replace('.fits','_base1.fits')
        file_new=first_look.baseline( file_in, file_out, 
                                      index_clean=index_rms, polyorder=1)
        first_look.peak_rms( file_new, index_rms=index_rms, 
                             index_peak=index_peak)

def write_skycoord_table(data, cube_ref, **kwargs):
    """
    Writes out a text file with flattened coordinates of the cube
    stacked with input array data. Additional arguments are passed
    to astropy's text writing function.

    TODO: add a useful `names` keyword?

    See astropy.io.ascii.write docstring for more info.

    Parameters
    ----------
    data : array-like structure of the same xy-grid as cube_ref.

    cube_ref : a cube file to get the coordinate grid from.

    """
    from astropy.table import Table
    from astropy.io import ascii
    from spectral_cube import SpectralCube

    cube = SpectralCube.read(cube_ref)

    flat_coords = [cube.spatial_coordinate_map[i].flatten() for i in [1,0]]
    # TODO: finish this up for multiple components
    #n_repeat = np.prod(np.array(data).shape)%np.prod(cube.shape[1:])+1

    table = Table(np.vstack(flat_coords +
        [np.array(xy_slice).flatten() for xy_slice in data]).T)

    ascii.write(table, **kwargs)

def moments(cube_fits, line_values, line_names, moment, save_file=False):
    """
    cube: str
        The datacube in fits format to open

    line_values: list of floats
        The wavelengths of the lines. !!! In general: if moment=0 the required wavelenth should be air, if moment=1 it should be vacuum !!!

    line_names: list of str
        The identifier of the lines

    moment: 0 or 1

    save_file: bool, optional
        Set to True if the result is to be saved as a fits file. Default is False.

    example:
        moment = moments('cube.fits', [4861.33, 6562.8], ['Hb', 'Ha'], moment=0)
    """


    print line_values, line_names
    cube=SpectralCube.read(cube_fits)

    for line,stri in zip(line_values,line_names):
        if moment==0:
            mom = cube.spectral_slab((line-3)*u.AA, (line+3)*u.AA).sum(axis=0)
            if save_file==True:
                mom.hdu.writeto(str(stri)+'_moment0.fits',clobber=True)
        if moment==1:    
            mom = cube.with_spectral_unit(u.km/u.s, rest_value=line*u.AA,velocity_convention='optical').spectral_slab(-300*u.km/u.s,300*u.km/u.s).moment1()
            if save_file==True:
                mom.hdu.writeto(str(stri)+'_moment1.fits',clobber=True)

    return mom

def S2_drawM33(vmin=40,vmax=80, deltaX=40, deltaV=6, deltadeltaX=10, deltadeltaV=1):
	"""Activates S2_draw with each of the .py file's subcube selections,
	   with the same args as S2_arrayM33.

	   Argument format: "(vmin=40,vmax=80, deltaX=40, deltaV=6, deltadeltaX=10,
	   deltadeltaV=1).

	   These MUST match the args/kwargs used in S2_arrayM33!"""

	galaxyname = 'M33'
	filename = 'm33.co21_iram_CLEANED'

	cube = SpectralCube.read(filename+".fits")

	pixelwidthDEG = cube.header['CDELT2']			# The width of each pixel, in degrees.
	distancePC = 840000.0			# The distance to the galaxy that M51's .fits file deals with, in parsecs.  (???) Is this number accurate, though?
	pixelwidthPC = pixelwidthDEG*np.pi/180.0*distancePC	# The width of each pixel, in pc.


	ymin = np.array([350,600,650,525,300,250])	# These are the minimum "y" values of the regions that we're dealing with.
	ymax = np.array([550,800,850,725,500,450])	# These are the corresponding maximum "y" values of these regions.
	xmin = np.array([500,100,400,288,200,550])	# These are the corresponding minimum "x" values of these regions.
	xmax = np.array([700,300,600,488,400,750])	# These are the corresponding maximum "x" values of these regions. (Example: The first region has ymin=350, ymax=550, xmin=500, xmax=700.)
	sets = np.ravel(ymin.shape)[0]		# This is the number of regions that we're dealing with.

	for i in range(0,sets):
		S2_draw(vmin,vmax,ymin[i],ymax[i],xmin[i],xmax[i],deltaX,deltaV,deltadeltaX,deltadeltaV,filename,galaxyname)

def cube_w11(region='IC348'):
    if region == 'IC348':
        OneOneFile = 'IC348mm/IC348mm-11_cvel_clean_rob05.fits'
        TwoTwoFile = 'IC348mm/IC348mm-11_cvel_clean_rob05.fits'
        vmin=7.4
        vmax=10.0
    elif region == 'IRAS03282':
        OneOneFile = 'IRAS03282/IRAS03282-11_cvel_clean_rob05.fits'
        TwoTwoFile = 'IRAS03282/IRAS03282-11_cvel_clean_rob05.fits'
        vmin=6.0
        vmax=8.5
    elif region == 'L1451mm':
        OneOneFile = 'L1451mm/L1451MM-11_cvel_clean_rob05.fits'
        TwoTwoFile = 'L1451mm/L1451MM-11_cvel_clean_rob05.fits'
        vmin=3.2
        vmax=4.9
    cube = SpectralCube.read(OneOneFile)
    vcube = cube.with_spectral_unit(u.km/u.s, rest_value=freq11, velocity_convention='radio')
    slab = vcube.spectral_slab( vmax*u.km/u.s, vmin*u.km/u.s)
    w11=slab.moment( order=0, axis=0)
    #beam = Beam.from_fits_header(fits.getheader(OneOneFile))
    # Next line is to solve bug in spectralcube:
    # it should be something like this in line 2234 of spectral_cube.py:
    # ```
    # if axis == 0 and self._meta['beam'] is not None:
    #     meta = { blabla, 'beam':self._meta['beam']}
    # else:
    #     meta = { blabla}
    w11._meta['beam'] = slab.beam
    w11.write(OneOneFile.replace('.fits','_w11.fits'), overwrite=True)

def writeplanes(save_name='/mnt/work/erosolow/GRS_13CO_all.fits'):
    spatial_template = fits.open('INTEG/COHRS_RELEASE1_FULL_INTEG.fit')
    spectral_template = SpectralCube.read('reprojected.fits')

    # Smoosh astrometry components together
    spatial_header = spatial_template[0].header
    spectral_header = spectral_template.header

    new_header = spatial_header.copy()
    new_header["NAXIS"] = 3
    for keyword in ['NAXIS3', 'CRVAL3', 'CDELT3','CRPIX3','CUNIT3']:
        new_header[keyword] = spectral_header[keyword]
    new_header['BMAJ'] = 14./3600
    new_header['BMIN'] = 14./3600
    new_header['BPA'] = 0.00
    
    if os.path.exists(save_name):
        raise Exception("The file name {} already "
                        "exists".format(save_name))

    # Open a file and start filling this with planes.
    output_fits = fits.StreamingHDU(save_name, new_header)
    # Again, set  up a common vel axis and spin out
    vel = np.linspace(-30, 160, 191)
    for v in vel:
        output_fits.write(fits.getdata(planesdir +
                                       'GRSPLANE_{0}'.format(v) +
                          '.fits'))

    output_fits.close()

def subcubes_from_ds9(cube, region_file='../nro_maps/SouthShells.reg', pad_factor=1., shape='exact'):
    """
    Extracts subcubes using the ds9 region file.
    
    Parameters
    ----------
    cube : SpectralCube, str
        The cube to be chopped. Must be type spectral_cube.SpectralCube or str filename.
    region_file : str
        Path to a ds9 region file.
    pad_factor : float, optional
        Expand the subcube around the region by this factor.
    shape : {'square', 'exact'}
        The shape of the subcube returned. 'square' returns the
        smallest square subcube that contains the region.
        'exact' returns only the pixels contained within the region.
    
    Returns
    -------
    subcubes: list of SpectralCube of SpectralCube
    """
    from spectral_cube import SpectralCube
    import pyregion

    try:
        #If cube is a str filename, read a SpectralCube.
        cube = SpectralCube.read(cube)
    except ValueError:
        pass

    if shape == 'square':
        import astropy.units as u
        subcube_list = []
        region_list = pyregion.open(region_file)
        for region in region_list:
            half_width = region.coord_list[2] * pad_factor * u.deg
            ra_center = region.coord_list[0] * u.deg
            dec_center = region.coord_list[1] * u.deg
            ra_range = [ra_center - half_width, ra_center + half_width]
            dec_range = [dec_center - half_width, dec_center + half_width]
            #print(ra_range, dec_range)
            subcube_list.append(cube.subcube(ra_range[1], ra_range[0], dec_range[0], dec_range[1]))
    if shape == 'exact':
        region_list = pyregion.open(region_file)
        subcube_list = []
        for region in region_list:
            
            if pad_factor != 1.:
                new_string = '{};{}({},{},{}")'.format(region.coord_format, region.name,
                                region.coord_list[0], region.coord_list[1],
                                region.coord_list[2]*3600.*pad_factor)
                region = pyregion.parse(new_string)[0]
                
            subcube_list.append(cube.subcube_from_ds9region(pyregion.ShapeList([region])))
    if len(subcube_list) == 1:
        return subcube_list[0]
    else:
        return subcube_list

def summary_plot(filelist):
    for thisfile in filelist:
        s = SpectralCube.read(thisfile)
        outfile = thisfile.replace('.fits','_summary.png')
        mom0 = s.moment0()
        f = aplpy.FITSFigure(mom0.hdu)
        f.show_colorscale()
        f.show_colorbar()
        f.save(outfile)

def cubegen(vmin,vmax,ymin,ymax,xmin,xmax, filename = "paws_norot", drawmap = False, mapname="3Dcube"):
	"""
	Returns a subcube of the specified dimensions from the .fits file.
	Also displays the subcube as it appears on the galaxy map if drawmap=True.


	Parameters:
	-----------
	vmin,...,xmax : int
		Parameters used in relevant xi map.
		WARNING: Selecting too large of a vmax-vmin will hugely increase
		processing time in later calculations.
	filename : str
		Name of the .paws data file.
		"paws_norot" for M51, "m33.co21_iram_CLEANED" for M33.
	drawmap : bool
		Enables or disables drawing the subcube Tmax map.
	galaxyname : str
		Name of the galaxy.
		'M51' for M51, 'M33' for M33.
	mapname : str
		Name of the saved image of the subcube's Tmax map, if
		drawmap==True.


	Returns:
	-----------
	subcube : spectral cube (?)
		The data inside the selected subcube.
"""

	cube = SpectralCube.read(filename+".fits")
	data = cube.filled_data[:]   # Pulls "cube"'s information (position, spectral info (?)) into a 3D Numpy array.
	yshape = data.shape[1]/2.0
	xshape = data.shape[2]/2.0

	pixelwidthDEG = cube.header['CDELT2']			# The width of each pixel, in degrees.
	if (filename =='m33.co21_iram_CLEANED') or (filename =='m33.co21_iram_CLEANED_smooth') or (filename =='m33.co21_iram_CLEANED_blank'):	# Checks if the galaxy's Header file contains its distance.
		distancePC = 840000.0				# The distance to the galaxy that M33's .fits file deals with, in parsecs. ONLY works on the CLEANED file!
	else:
		distancePC = cube.header['DIST']		# The distance to the galaxy that M51's .fits file deals with, in parsecs.  (???) Is this number accurate, though?
	pixelwidthPC = pixelwidthDEG*np.pi/180.0*distancePC	# The width of each pixel, in pc.

	subcube = cube[vmin:vmax,ymin:ymax,xmin:xmax]
	if drawmap == True:
		plt.figure(1)
		plt.imshow(np.nanmax(data[vmin:vmax,ymin:ymax,xmin:xmax].value,axis=0), extent=[(xmin-xshape)*pixelwidthPC,(xmax-xshape)*pixelwidthPC, \
			   (ymin-yshape)*pixelwidthPC,(ymax-yshape)*pixelwidthPC], origin='lower')
		fig = matplotlib.pyplot.gcf()
		#fig.set_size_inches(5, 5)	# Enlarges the image so as to prevent squishing.
		plt.xlabel('Distance from Centre in x-direction (pc)')
		plt.ylabel('Distance from Centre in y-direction (pc)')

		plt.savefig('galaxy_'+mapname+'.png')
		plt.clf()			# Clears the image after saving.

	return subcube

def FirstLook_NGC1333():
    print("Now NH3(1,1)")
    a_rms = [  0, 158, 315, 428, 530, 693, 751]
    b_rms = [ 60, 230, 327, 438, 604, 735, 760]
    index_rms=first_look.create_index( a_rms, b_rms)
    index_peak=np.arange(326,430)
    file_in='NGC1333/NGC1333_NH3_11.fits'
    # 1st order polynomial
    file_out=file_in.replace('.fits','_base1.fits')
    file_new=first_look.baseline( file_in, file_out, index_clean=index_rms, polyorder=1)
    first_look.peak_rms( file_new, index_rms=index_rms, index_peak=index_peak)
    print("Now NH3(2,2)")
    a_rms = [  0, 190, 360, 600]
    b_rms = [70, 300, 470, 640]
    index_rms=first_look.create_index( a_rms, b_rms)
    index_peak=np.arange(380,520)
    # file_in='NGC1333/NGC1333_NH3_22.fits'
    # # 1st order polynomial
    # file_out=file_in.replace('.fits','_base1.fits')
    # file_new=first_look.baseline( file_in, file_out, index_clean=index_rms, polyorder=1)
    # first_look.peak_rms( file_new, index_rms=index_rms, index_peak=index_peak)
    ## 2nd order polynomial
    #file_out=file_in.replace('.fits','_base2.fits')
    #file_new=first_look.baseline( file_in, file_out, index_clean=index_rms, polyorder=2)
    #first_look.peak_rms( file_new, index_rms=index_rms, index_peak=index_peak)
    #
    # print("Now NH3(3,3)")
    # a_rms = [ 10, 190, 420]
    # b_rms = [70, 360, 500]
    # index_rms=first_look.create_index( a_rms, b_rms)
    # index_peak=np.arange(410,540)
    # file_in='NGC1333/NGC1333_NH3_33.fits'
    # 1st order polynomial
    # file_out=file_in.replace('.fits','_base1.fits')
    # file_new=first_look.baseline( file_in, file_out, index_clean=index_rms, polyorder=1)
    # first_look.peak_rms( file_new, index_rms=index_rms, index_peak=index_peak)
    linelist = ['NH3_22','NH3_33','C2S','HC5N','HC7N_21_20','HC7N_22_21']
    vsys = 8.5*u.km/u.s
    throw = 8*u.km/u.s
    for line in linelist:
        file_in = 'NGC1333/NGC1333_{0}.fits'.format(line)
        s = SpectralCube.read(file_in)
        s = s.with_spectral_unit(u.km/u.s,velocity_convention='radio')
        a_rms = [s.closest_spectral_channel(vsys+2*throw),
                 s.closest_spectral_channel(vsys-throw)]
        b_rms = [s.closest_spectral_channel(vsys+throw),
                 s.closest_spectral_channel(vsys-2*throw)]
        index_peak = np.arange(s.closest_spectral_channel(vsys+3*u.km/u.s),
                              s.closest_spectral_channel(vsys-3*u.km/u.s))
        index_rms=first_look.create_index( a_rms, b_rms)

        file_out=file_in.replace('.fits','_base1.fits')
        file_new=first_look.baseline( file_in, file_out, 
                                      index_clean=index_rms, polyorder=1)
        first_look.peak_rms( file_new, index_rms=index_rms, 
                             index_peak=index_peak)

def dendropix(fileprefix='SgrB2_b3_12M.HC3N'):
    cube = SpectralCube.read(dpath('{0}.image.pbcor.contsub.fits'.format(fileprefix))).minimal_subcube()
    noise = cube.spectral_slab(-200*u.km/u.s, -100*u.km/u.s).std(axis=0)
    keep_mask = cube.max(axis=0) > noise

    tblfile = tpath('{0}.dendrotable.ecsv'.format(fileprefix))
    if os.path.exists(tblfile):
        table = Table.read(tblfile, format='ascii.ecsv')
    else:
        table = Table([Column(name='xpix'),
                       Column(name='ypix'),
                       Column(name='zpix'),
                       Column(name='lon'),
                       Column(name='lat'),
                       Column(name='velo'),
                       Column(name='peakval'),])

    xpyp_done = set(zip(table['xpix'], table['ypix']))
    all_keepers = zip(*np.where(keep_mask))
    xpyp = [x for x in all_keepers if x not in xpyp_done]
    print(len(xpyp), len(all_keepers), len(xpyp_done))
    if len(xpyp_done) > 0:
        assert len(xpyp) < len(all_keepers)


    for ii,(ypix,xpix) in enumerate(ProgressBar(xpyp)):

        data = cube[:,ypix,xpix].value

        error = noise[ypix,xpix].value
        # alternative:
        #error = stats.sigma_clipped_stats(data)[2]

        D = astrodendro.Dendrogram.compute(data, min_value=0,
                                           min_delta=2*error, min_npix=7,
                                           is_independent=astrodendro.pruning.min_peak(5*error))
        if not D.leaves:
            table.add_row([xpix,ypix,]+[np.nan]*5)
            del D
            continue

        #peaks = [S.get_peak()[0][0] for S in D]
        #peak_vals = [S.get_peak()[1] for S in D]
        #peaks = [cube.spectral_axis[S.get_peak()[0][0]].to(u.km/u.s).value for S in D]

        for S in D:
            (peak_pix,),peak_val = S.get_peak()
            velo,lat,lon = cube.world[peak_pix, ypix, xpix]
            table.add_row([xpix,ypix,peak_pix,lon,lat,velo,peak_val])

        if ii % 100 == 0:
            table.write(tblfile, format='ascii.ecsv')

        del D
        del S
        del data