LTE Line Model example: Line identification

import numpy as np

# retrieve the data
import requests

# 10 MB file: stream it, but it should be able to load in memory
response = requests.get('', stream=True)

with open('temp_cube.fits', 'wb') as fh:

from spectral_cube import SpectralCube 
from astropy import units as u
# there is some junk at low frequency
cube ='temp_cube.fits').spectral_slab(218.14*u.GHz, 218.54*u.GHz)

# this isn't really right b/c of beam issues, but it's Good Enough
meanspec = cube.mean(axis=(1,2))
# have to hack this, which is _awful_ and we need to fix it
meanspec_K = meanspec.value * cube.beam.jtok(cube.spectral_axis)

import pyspeckit

sp = pyspeckit.Spectrum(xarr=cube.spectral_axis, data=meanspec_K)

# The LTE model doesn't include an explicit filling factor, so we impose one
fillingfactor = 0.05

# there's continuum in our spectrum, but it's also not modeled
offset = np.nanpercentile(, 20)

# look at a couple species
# Be cautious! Not all line catalogs have all of these species, some species
# can result in multiple matches from splatalogue, and definitely not all
# species have the same column and excitation
species_list = ('CH3OH', 'HCCCN', 'H2CO', 'CH3OCHO', 'CH3OCH3', 'C2H5CN',)
# 'C2H3CN' is pretty definitely not here at the constant assumed column

# make a set of subplots:
import pylab as pl
fig, axes = pl.subplots(len(species_list)+1, 1, figsize=(18,14), num=1)

from pyspeckit.spectrum.models import lte_molecule
mods = []
for species, axis in zip(species_list, axes):
    # search for lines within +/- 0.5 GHz of the min/max (to account for velocity offset, generously)
    freqs, aij, deg, EU, partfunc = lte_molecule.get_molecular_parameters(species,

    mod = lte_molecule.generate_model(xarr=sp.xarr, vcen=50*,
                                      width=3*, tex=100*u.K,
                                      column=1e17***-2, freqs=freqs,
                                      aij=aij, deg=deg, EU=EU,

    sp.plotter(axis=axis, figure=fig, clear=False)
    sp.plotter.axis.plot(sp.xarr, mod*fillingfactor+offset, label=species, zorder=-1)

    if axis != axes[-1]:
    if axis != axes[len(species_list) // 2]:

axis = axes[-1]
sp.plotter(axis=axis, figure=fig, clear=False)
sp.plotter.axis.plot(sp.xarr, np.nansum(mods,axis=0)*fillingfactor+offset, label='Sum', zorder=-1)
W51 spectrum with line IDs

W51 spectrum with LTE lines at T_X = 100, N=10^17

[edit this page on github]

[edit this page on bitbucket]