用spectral_cube抽取cube文件上某个点或者面源的光谱

from astropy.io import fits
from spectral_cube import SpectralCube as spec
import numpy as np
import matplotlib.pyplot as plt
import os
import sys
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
from astropy.table import Table
import warnings
warnings.filterwarnings('ignore')
T = Table.read('atoms_B3_12M_allsour_cont.tbl',format='ascii')
#把需要的位置信息列成表格，表格格式如下


def seek_file(sourcename):
    files = []
    for filename in os.listdir(os.getcwd()):
        fp = os.path.join(filename)
        spw35 = sourcename+'_sci.spw35.cube.I.pbcor.fits'
        spw37 = sourcename+'_sci.spw37.cube.I.pbcor.fits'
        spw39 = sourcename+'_sci.spw39.cube.I.pbcor.fits'
        if os.path.isfile(fp) and spw35 in filename:
            files.append(fp)
        if os.path.isfile(fp) and spw37 in filename:
            files.append(fp)
        if os.path.isfile(fp) and spw39 in filename:
            files.append(fp)
    if len(files)==3:
        files.sort()
        files = files[1:3]
    return files

def extra(sourcename,files,ra,dec):
    import astropy.units as u
    cube = spec.read(files[0])
#    cube.allow_huge_operations=True
#    cube = cube.to(u.K)
    hdu = fits.open(files[0])
    hdr = hdu[0].header
    delta1 = hdr['CDELT1']*3600
    delta2 = hdr['CDELT2']*3600
    offset1 = int(hdr['CRPIX1']-(hdr['crval1']-ra)*3600/delta1*np.cos(dec/180*np.pi))
    offset2 = int(hdr['CRPIX2']+(hdr['crval2']-dec)*3600/delta1)

    subcube = cube[:,offset2-4:offset2+4,offset1-4:offset1+4]
    spectrum = subcube.mean(axis=(1,2))
#    bmaj = spetrum.header['BMAJ']*3600*u.arcsec
#    bmin = spetrum.header['BMIN']*3600*u.arcsec
#    fwhm_to_sigma = 1./(8*np.log(2))**0.5
#    beam_area = 2.*np.pi*(bmaj*bmin*fwhm_to_sigma**2)
#    freq = 98*u.GHz
#    equiv = u.brightness_temperature(freq)
    spectrum = spectrum.to(u.K)
    coord = SkyCoord(ra,dec,unit='deg',frame='fk5')
    hh,hm,hs = coord.ra.hms
    dd,dm,ds = coord.dec.dms
    dm,ds = abs(dm),abs(ds)
    spectrum.write('spec/%s_%ih%im%.4fs_%id%im%.4fs_98G.fits' %(sourcename,hh,hm,hs,dd,dm,ds))

    cube = spec.read(files[1])
    hdu = fits.open(files[1])
    hdr = hdu[0].header
    delta1 = hdr['CDELT1']*3600
    delta2 = hdr['CDELT2']*3600
    offset1 = int(hdr['CRPIX1']-(hdr['crval1']-ra)*3600/delta1*np.cos(dec/180*np.pi))
    offset2 = int(hdr['CRPIX2']+(hdr['crval2']-dec)*3600/delta1)

    subcube = cube[:,offset2-4:offset2+4,offset1-4:offset1+4]
    spectrum = subcube.mean(axis=(1,2))
    spectrum = spectrum.to(u.K)
    coord = SkyCoord(ra,dec,unit='deg',frame='fk5')
    hh,hm,hs = coord.ra.hms
    dd,dm,ds = coord.dec.dms
    dm,ds = abs(dm),abs(ds)
    spectrum.write('spec/%s_%ih%im%.4fs_%id%im%.4fs_100G.fits' %(sourcename,hh,hm,hs,dd,dm,ds))

for i in range(len(T['sname'])):
    files = seek_file(T['sname'][i])
    print(files)
    if len(files)==2:
        extra(T['sname'][i],files,T['ra'][i],T['dec'][i])