Preparation of Data

Before starting the line fit we need to prepare the data. As an example, if different lines need to be combined later in the analysis, then my suggestion is to first prepare all the cubes to a common beam and grid.

Convolving to a common beam

The following code shows how to convolve a cube to a common beam. It is implemented to work witht a list of files to be convolved. The code will convolve the cubes to the beam of the first cube in the list.

 1import radio_beam
 2import astropy.units as u
 3from astropy.io import fits
 4from spectral_cube import SpectralCube
 5
 6# takes 2 files, and find the common beam
 7# and convolve both to that beam
 8# then write out the new files
 9
10# files to be convolved
11fits_files = ['cube1.fits', 'cube2.fits']
12# suffix for the new files
13suffix = '_convolved'
14
15# load the parameters of the beams
16bmaj_list, bmin_list, bpa_list = [], [], []
17for file_i in fits_files:
18    hd = fits.getheader(file_i)
19    bmaj_list.append(hd['BMAJ'])
20    bmin_list.append(hd['BMIN'])
21    bpa_list.append(hd['BPA'])
22
23# create the Beam objects
24my_beams = radio_beam.Beams(
25    major=bmaj_list * u.deg, minor=bmin_list * u.deg, pa=bpa_list * u.deg)
26common_beam = my_beams.common_beam()
27
28# write out the smoothed cubes
29for file_in in fits_files:
30    cube = SpectralCube.read(file_in)
31    new_cube = cube.convolve_to(common_beam)
32    new_cube.write(file_in.replace(
33        '.fits', f'{suffix}.fits'), overwrite=True)

File regridding

The regridding of cubes to a matched grid can be done by leveraging the options already available in the spectral-cube package. The following code shows how to regrid a cube to a common grid.

 1from astropy.io import fits
 2from spectral_cube import SpectralCube
 3
 4# files to be regridded
 5fits_files = ['cube1.fits', 'cube2.fits']
 6# suffix for the new files
 7suffix = '_regridded'
 8# template file
 9template_file = 'template.fits'
10
11# load the template cube
12# check if the template file has the necessary header keys
13# there might be more keys that need to be checked
14template_hd = fits.getheader(template_file)
15list_key = ['NAXIS1', 'NAXIS2', 'CRVAL1', 'CRVAL2', 'CDELT1', 'CDELT2',
16            'CRPIX1', 'CRPIX2', 'CTYPE1', 'CTYPE2', 'CUNIT1', 'CUNIT2',]
17for key in list_key:
18    if key not in template_hd:
19        raise KeyError(f'Header key {key} not found in the template file.')
20# load all cubes and regrid them
21for f in fits_files:
22    cube = SpectralCube.read(f)
23    target_hd = cube.header
24    for key in list_key:
25        target_hd[key] = template_hd[key]
26    regridded_cube = cube.reproject(target_hd)
27    regridded_cube.write(
28        f.replace('.fits', f'{suffix}.fits'), overwrite=True)

Unit conversion and primary beam correction

The convertion of interferometric data into Kelvin units and including the primary beam conversion can be performed with the following code.

 1from spectral_cube import SpectralCube
 2import numpy as np
 3import astropy.units as u
 4
 5# file without the primary beam corrected
 6# then primary beam response, and primary beam corrected one
 7file_in = 'cube_no_PBcorr.fits'
 8file_in_PB = 'cube_PB_response.fits'
 9file_in_K = 'cube_ready.fits'
10
11cube = SpectralCube.read(file_in)
12PB = fits.getdata(file_in_PB)
13kcube = cube.to(u.K) / np.squeeze(PB)
14kcube.write(file_in_K, overwrite=True)