Python SpectralCube.read方法代码示例

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import read [as 别名]
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

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import read [as 别名]
    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)

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import read [as 别名]
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

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import read [as 别名]
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)

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import read [as 别名]
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)

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import read [as 别名]
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)

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import read [as 别名]
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)

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import read [as 别名]
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)

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import read [as 别名]
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)

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import read [as 别名]
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

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import read [as 别名]
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()

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import read [as 别名]
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

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import read [as 别名]
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

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import read [as 别名]
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)

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import read [as 别名]
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)