Linetools

linetools is an in-development package for the analysis of 1-D spectra, with the aim to become an Astropy affiliated package. Its core developers work primarily on UV/optical/IR absorption line research, so most of the functionality is aimed at the identification and analysis of absorption lines. The eventual goal is to provide a set of tools useful for both absorption and emission lines.

Note

linetools is still under active development. While the developers strive to maintain compatibility in new releases, there may be backwards-incompatible changes in future versions.

Getting Started

Installation

Dependencies

Linetools depends on these packages:

  • python versions 2.7, or 3.3 or later
  • numpy version 1.9 or later
  • astropy version 1.0 or later
  • scipy version 0.16 or later
  • matplotlib version 1.4 or later
  • PyQT4 version 4 (for GUIs)
  • h5py version 2.6 (for data I/O)

We strongly recommend that you use Anaconda to install them. With Anaconda you can check for the presence and versions of the dependencies with:

conda list "^python|numpy|astropy|scipy|matplotlib|PyQT|h5py"

If you’re missing any, install them with (for example):

conda install astropy scipy matplotlib PyQT h5py

If their versions are too old, update them with (for example):

conda update astropy

If you aren’t using Anaconda then all of the dependencies can also be installed with pip.

Installing Linetools

If you plan to play around with the code and possibly contribute changes, then follow the instructions in the section below, Installing Linetools from Source. Otherwise simply use:

pip install linetools

and you’re done!

If you wish to have full functionality of the GUIs and are using MacOSX, then you probably need to change your backend from macosx to TkAgg in the matplotlibrc file.

Installing Linetools from Source

I just want to play with the code

Install the development version like this:

git clone https://github.com/linetools/linetools.git
cd linetools
python setup.py develop

Now you can easily make tweaks to the code, which are immediately applied to your installed version (you’ll have to reload the relevant modules to see those changes in an existing Python session, though).

I want to make a code contribution to linetools

Fantastic! In that case, follow the Astropy developer guidelines, replacing every instance of astropy in those instructions with linetools. This will install a ‘fork’ of linetools that you can use to submit code changes to the main repository.

Running Tests

To test your installation, run:

python -c 'import linetools; linetools.test()'

The tests take a couple of minutes to finish. If you notice any failures, we’d love you to report them on the linetools issue tracker.

Before Launching GUIs

If you are a Mac user, we highly recommend that you set your matplotlib backend from MacOSX to TkAgg (or another option, see backends).

Building Documentation

Only do this if you’re a developer! If you want build the documentation, you also need to install Sphinx (version 1.3+):

conda install sphinx

If you’d like to generate inheritance diagrams in the docs then you also need to install graphviz (MacOSX, Ubuntu), but this isn’t required. Once sphinx is installed, change to the /docs directory under the source directory and run:

make html

The documentation should now be in _build/html.

Changelog

0.2 (unreleased)

Updates
  • LineList.available_transitions() no longer has key argument n_max
  • LineList: extra attributes for transitions added (ion_name, log(w*f), abundance, ion_correction, rel_strength)
  • ASCII tables with no header are required to be 4 columns or less for io.readspec to work
  • Modify XSpectrum1D to use masked numpy arrays
  • Enable hdf5 I/O [requires h5py]
  • Added .header property to XSpectrum1D (reads from .meta)
  • Added XSpectrum1D.write, a generic write wrapper
  • Added xabssysgui GUI
  • Added new linelists (e.g. Galaxy)
  • Added EmissLine child to SpectralLine
  • Added LineLimits class
  • Added SolarAbund class
  • Added lt_radec and lt_line scripts
Bug fixes
  • Fix XSpectrum1D.from_tuple to allow an astropy table column to specify wavelengths and fluxes.

0.1 (2016-01-31)

First public release.

Core classes

AbsComponent Class

Notebooks

Examples for the AbsComponent Class (v0.3)

Download this notebook.

%matplotlib inline

# suppress warnings for these examples
import warnings
warnings.filterwarnings('ignore')

# import
try:
    import seaborn as sns; sns.set_style("white")
except:
    pass

import astropy.units as u
from linetools.spectralline import AbsLine
from linetools.isgm import utils as ltiu
from linetools.analysis import absline as laa
from linetools.spectra import io as lsio
from linetools.isgm.abscomponent import AbsComponent

import imp
lt_path = imp.find_module('linetools')[1]
Instantiate
Standard 
abscomp = AbsComponent((10.0*u.deg, 45*u.deg), (14,2), 1.0, [-300,300]*u.km/u.s)
abscomp

<AbsComponent: 00:40:00 +45:00:00, Name=SiII_z1.00000, Zion=(14,2), Ej=0 1 / cm, z=1, vlim=-300 km / s,300 km / s>
From AbsLines
From one line 
lya = AbsLine(1215.670*u.AA)
lya.analy['vlim'] = [-300.,300.]*u.km/u.s
lya.attrib['z'] = 2.92939

abscomp = AbsComponent.from_abslines([lya])
print(abscomp)
abscomp._abslines

<AbsComponent: 00:00:00 +00:00:00, Name=HI_z2.92939, Zion=(1,1), Ej=0 1 / cm, z=2.92939, vlim=-300 km / s,300 km / s>

[<AbsLine: HI 1215, wrest=1215.6700 Angstrom>]
From multiple 
lyb = AbsLine(1025.7222*u.AA)
lyb.analy['vlim'] = [-300.,300.]*u.km/u.s
lyb.attrib['z'] = lya.attrib['z']

abscomp = AbsComponent.from_abslines([lya,lyb])
print(abscomp)
abscomp._abslines

<AbsComponent: 00:00:00 +00:00:00, Name=HI_z2.92939, Zion=(1,1), Ej=0 1 / cm, z=2.92939, vlim=-300 km / s,300 km / s>

[<AbsLine: HI 1215, wrest=1215.6700 Angstrom>,
 <AbsLine: HI 1025, wrest=1025.7222 Angstrom>]
Methods
Generate a Component Table 
lya.attrib['logN'] = 14.1
lya.attrib['sig_logN'] = 0.15
lya.attrib['flag_N'] = 1
laa.linear_clm(lya.attrib)
lyb.attrib['logN'] = 14.15
lyb.attrib['sig_logN'] = 0.19
lyb.attrib['flag_N'] = 1
laa.linear_clm(lyb.attrib)

(<Quantity 141253754462275.53 1 / cm2>, <Quantity 61797269977312.6 1 / cm2>)

abscomp = AbsComponent.from_abslines([lya,lyb])
comp_tbl = abscomp.build_table()
comp_tbl
<QTable length=2>
wrestzflag_NlogNsig_logN
Angstrom
float64float64int64float64float64
1215.672.92939114.10.15
1025.72222.92939114.150.19
Synthesize multiple components 
SiIItrans = ['SiII 1260', 'SiII 1304', 'SiII 1526']
SiIIlines = []
for trans in SiIItrans:
    iline = AbsLine(trans)
    iline.attrib['logN'] = 12.8 + np.random.rand()
    iline.attrib['sig_logN'] = 0.15
    iline.attrib['flag_N'] = 1
    iline.attrib['z'] = 2.92939
    iline.analy['vlim'] = [-300.,50.]*u.km/u.s
    _,_ = laa.linear_clm(iline.attrib)
    SiIIlines.append(iline)
SiIIcomp = AbsComponent.from_abslines(SiIIlines)
SiIIcomp.synthesize_colm()

SiIIlines2 = []
for trans in SiIItrans:
    iline = AbsLine(trans)
    iline.attrib['logN'] = 13.3 + np.random.rand()
    iline.attrib['sig_logN'] = 0.15
    iline.attrib['flag_N'] = 1
    iline.attrib['z'] = 2.92939
    iline.analy['vlim'] = [50.,300.]*u.km/u.s
    _,_ = laa.linear_clm(iline.attrib)
    SiIIlines2.append(iline)
SiIIcomp2 = AbsComponent.from_abslines(SiIIlines2)
SiIIcomp2.synthesize_colm()

abscomp.synthesize_colm()
[abscomp,SiIIcomp,SiIIcomp2]

[<AbsComponent: 00:00:00 +00:00:00, Name=HI_z2.92939, Zion=(1,1), Ej=0 1 / cm, z=2.92939, vlim=-300 km / s,300 km / s, logN=14.1172, sig_N=0.117912, flag_N=1>,
 <AbsComponent: 00:00:00 +00:00:00, Name=SiII_z2.92939, Zion=(14,2), Ej=0 1 / cm, z=2.92939, vlim=-300 km / s,50 km / s, logN=12.9226, sig_N=0.112727, flag_N=1>,
 <AbsComponent: 00:00:00 +00:00:00, Name=SiII_z2.92939, Zion=(14,2), Ej=0 1 / cm, z=2.92939, vlim=50 km / s,300 km / s, logN=13.8523, sig_N=0.0897197, flag_N=1>]

synth_SiII = ltiu.synthesize_components([SiIIcomp,SiIIcomp2])
synth_SiII

<AbsComponent: 00:00:00 +00:00:00, Name=SiII_z2.92939, Zion=(14,2), Ej=0 1 / cm, z=2.92939, vlim=-300 km / s,300 km / s, logN=13.9006, sig_N=0.0811523, flag_N=1>
Generate multiple components from abslines 
comps = ltiu.build_components_from_abslines([lya,lyb,SiIIlines[0],SiIIlines[1]])
comps

[<AbsComponent: 00:00:00 +00:00:00, Name=HI_z2.92939, Zion=(1,1), Ej=0 1 / cm, z=2.92939, vlim=-300 km / s,300 km / s>,
 <AbsComponent: 00:00:00 +00:00:00, Name=SiII_z2.92939, Zion=(14,2), Ej=0 1 / cm, z=2.92939, vlim=-300 km / s,50 km / s>]
Generate an Ion Table 
tbl = ltiu.iontable_from_components([abscomp,SiIIcomp,SiIIcomp2])
tbl
<QTable length=2>
ZionAEjzvminvmaxflag_NlogNsig_logN
km / skm / s
int64int64int64float64float64float64float64int64float64float64
1100.02.92939-300.0300.0114.11720248170.117911610801
14200.02.92939-300.0300.0113.90061577330.0811522506077
Stack plot
Load a spectrum 
xspec = lsio.readspec(lt_path+'/spectra/tests/files/UM184_nF.fits')
lya.analy['spec'] = xspec
lyb.analy['spec'] = xspec
Show 
abscomp = AbsComponent.from_abslines([lya,lyb])
abscomp.stack_plot()
_images/AbsComponent_examples_30_0.png
Column Densities with AbsComponent

Download this notebook.

%matplotlib inline

# suppress warnings for these examples
import warnings
warnings.filterwarnings('ignore')

# imports
try:
    import seaborn as sns; sns.set(context="notebook",font_scale=2)
except:
    pass

from scipy import integrate
import astropy.units as u

from linetools.isgm import abscomponent as lt_abscomp
from linetools.spectralline import AbsLine
from linetools.spectra.xspectrum1d import XSpectrum1D
#
import imp
lt_path = imp.find_module('linetools')[1]
Read Spec
xspec = XSpectrum1D.from_file(lt_path+'/spectra/tests/files/UM184_nF.fits')
Generate a few AbsLines
SiIItrans = ['SiII 1260', 'SiII 1304', 'SiII 1526', 'SiII 1808']

abslines = []
for trans in SiIItrans:
    iline = AbsLine(trans)
    iline.attrib['z'] = 2.92939
    iline.analy['vlim'] = [-250.,80.]*u.km/u.s
    iline.analy['spec'] = xspec
    abslines.append(iline)
#
abslines

[<AbsLine: SiII 1260, wrest=1260.4221 Angstrom>,
 <AbsLine: SiII 1304, wrest=1304.3702 Angstrom>,
 <AbsLine: SiII 1526, wrest=1526.7070 Angstrom>,
 <AbsLine: SiII 1808, wrest=1808.0129 Angstrom>]
Generate the Component
abscomp = lt_abscomp.AbsComponent.from_abslines(abslines)

try:
    sns.set(context="notebook",font_scale=2)
except:
    pass
abscomp.stack_plot()
_images/AbsComponent_ColumnDensities_10_0.png
Synthesize/Measure AODM Column Densities
abscomp.synthesize_colm(redo_aodm=True)

abscomp.logN

13.594445560856554

for iline in abscomp._abslines:
    print(iline.wrest, iline.attrib['flag_N'], iline.attrib['logN'], iline.attrib['sig_logN'])

1260.4221 Angstrom 1 13.5883729709 0.0150745701489
1304.3702 Angstrom 1 13.7708705955 0.0862006463782
1526.707 Angstrom 1 13.6707360009 0.0640855113383
1808.0129 Angstrom 3 0.0 0.50976387151
Apparent Column Density Plot
abscomp.plot_Na()
_images/AbsComponent_ColumnDensities_18_0.png
COG
\(F(\tau_0)\)

Definition \(F(\tau_0) = \int_0^\infty dx \, [1- \rm e^{-\tau_0 \rm e^{-x^2}}]\)

def ftau_intgrnd(x,tau0=0.1):
    return 1 - np.exp(-tau0 * np.exp(-x**2))

neval = 10000
lgt = np.linspace(-3, 9, neval)
all_tau0 = 10.**lgt
Ftau = np.zeros(neval)
for jj,tau0 in enumerate(all_tau0):
    Ftau[jj], ferr = integrate.quad(ftau_intgrnd, 0, np.inf, args=(tau0,))

# Damped limit (not accurate enough)
damp_lgt = np.linspace(6, 10, 100)
damp_tau0 = 10.**damp_lgt
damp_Ftau = np.sqrt(np.log(damp_tau0))

import matplotlib.pyplot as plt
plt.plot(lgt, Ftau, damp_lgt, 1.015*damp_Ftau)

[<matplotlib.lines.Line2D at 0x10c48b5c0>,
 <matplotlib.lines.Line2D at 0x10c464e10>]
_images/AbsComponent_ColumnDensities_25_1.png
Perform and Plot 
abscomp = lt_abscomp.AbsComponent.from_abslines(abslines)
COG_dict = abscomp.cog(redo_EW=True, show_plot=True)
_images/AbsComponent_ColumnDensities_27_0.png 
# Output
COG_dict

{'EW': <Quantity [ 0.43129915, 0.06810455, 0.11137664,-0.01950807] Angstrom>,
 'b': <Quantity 49.22868767597288 km / s>,
 'f': array([ 1.18   ,  0.0863 ,  0.127  ,  0.00208]),
 'logN': 13.693355878125537,
 'parm': <single_cog_model(logN=13.693355878125537, b=49.22868767597288)>,
 'redEW': array([  3.42186280e-04,   5.22125891e-05,   7.29522068e-05,
         -1.07897867e-05]),
 'sig_EW': <Quantity [ 0.0129661 , 0.01440996, 0.01686854, 0.02102034] Angstrom>,
 'sig_b': <Quantity 6.356381185059458 km / s>,
 'sig_logN': 0.054323725737309987,
 'wrest': <Quantity [ 1260.4221, 1304.3702, 1526.707 , 1808.0129] Angstrom>}

Overview

This Class is designed to organize and analyze a set of absorption lines.

By definition, an AbsComponent is a unique collection of absorption lines specified by:

Property Variable Type Description
RA, DEC radec tuple or coord RA,DEC in deg or astropy.coordinate
Z, ion Zion tuple Atomic Number (Z) and ionization state (ion)
Redshift z float absorption redshift
Velocity limits vlim Quantity array -/+ velocity limits of the component
Energy level Ej Quantity Energy of the level relative to ground
Isotope A int Atomic Mass number (optional)

Instantiation

The AbsComponent Class may be instantiated in a few ways. The default sets the properties listed above:

abscomp = AbsComponent((10.0*u.deg, 45*u.deg), (14,2), 1.0, [-300,300]*u.km/u.s)

More commonly, one will instantiate with one AbsLine object:

lya = AbsLine(1215.670*u.AA, z=2.92939)
lya.limits.set([-300.,300.]*u.km/u.s)  # vlim
abscomp1 = AbsComponent.from_abslines([lya])

or multiple:

lyb = AbsLine(1025.7222*u.AA, z=lya.attrib['z'])
lyb.limits.set([-300.,300.]*u.km/u.s)  # vlim
abscomp = AbsComponent.from_abslines([lya,lyb])

One may also instantiate from a dict, usually read from the hard-drive:

abscomp = AbsComponent.from_dict(idict)

Inspecting

Here are a few simple methods to explore/inspect the class.

Generate a QTable

If the class contains one or more AbsLines, you may generate a QTable from their attributes and data:

comp_tbl = abscomp.build_table()
Show a Stack Plot

If the AbsLine have spectra attached to them (in attrib[‘spec’]), a stack plot (aka velocity plot) is generated with:

abscomp.stack_plot()
Apparent Column Densitities

Show a plot of the apparent column density profiles, \(N_a\):

abscomp.plot_Na()

Analysis

Here are some methods related to analysis.

Synthesize Columns

If one inputs a set of AbsLine(s) with column density measurements, the synthesize_colm method collates these. Positive, unsaturated detections are combined in a weighted mean whereas limits are compared and the most sensitive one is adopted.:

abscomp.synthesize_colm()

Here is the set of rules:

  1. If all measurements are upper limits, take the lowest value and flag as an upper limit (flgN=3).
  2. If all measurements are a mix of upper and lower limits, take the highest lower limit and flag as a lower limit (flgN=2).
  3. If one or more measurements are a proper detection, take the weighted mean of these and flag as a detection (flgN=1).
Curve of Growth

A standard, single-component curve-of-growth (COG) analysis may be performed on the set of AbsLines:

COG_dict = abscomp.cog(show_plot=True)

The output dict includes:

Key Type Description
EW Quantity array Input equivalent widths
sigEW Quantity array Input error in equivalent widths
f ndarray Input f-values
wrest Quantity array Input rest wavelengths
logN float Output fitted column density (log10)
sig_logN float Output error in fitted logN
b Quantity Output b-value (km/s)
sig_b Quantity Output error in b-value (km/s)

Misc

I/O

One may generate a dict of the key properties of the AbsSystem with the to_dict() method:

cdict = component.to_dict()
Synthesize Components

This method combines a list of two or more components into a new one. It checks first for consistent RA/DEC, Zion, and Ej. It does not place any constraints on z and vlim. The column density of the new component is the sum of the input ones (with rules for limits). And the redshift and vlim are set to encompass the velocities of the input components.:

from linetools.isgm import utils as ltiu
synth_SiII = ltiu.synthesize_components([SiIIcomp1,SiIIcomp2])

See the Examples for the AbsComponent Class (v0.3) notebook for a complete example.

Generate Multiple Components

This method generates multiple components from a list of AbsLines.:

comps = ltiu.build_components_from_abslines([lya,lyb,SiIIlines[0],SiIIlines[1]])

AbsSystem Class

Notebooks

Examples for AbsSystem Class (v1.1)

Download this notebook.

# suppress warnings for these examples
import warnings
warnings.filterwarnings('ignore')

# imports
import imp
from astropy.coordinates import SkyCoord
import astropy.units as u

from linetools.isgm import abssystem as lt_absys
from linetools.spectralline import AbsLine
from linetools.isgm.abscomponent import AbsComponent
Simple instantiation
Standard init 
radec = SkyCoord(ra=123.1143*u.deg, dec=-12.4321*u.deg)
gensys = lt_absys.GenericAbsSystem(radec, 1.244, [-500,500]*u.km/u.s, NHI=16.)
gensys

<GenericAbsSystem: name=Foo type=Generic, 08:12:27.432 -12:25:55.56, z=1.244, NHI=16>
From components
One component 
# HI Lya, Lyb
lya = AbsLine(1215.670*u.AA)
lya.analy['vlim'] = [-300.,300.]*u.km/u.s
lya.attrib['z'] = 2.92939
lyb = AbsLine(1025.7222*u.AA)
lyb.analy['vlim'] = [-300.,300.]*u.km/u.s
lyb.attrib['z'] = lya.attrib['z']
abscomp = AbsComponent.from_abslines([lya,lyb])
abscomp.coord = radec

linetools.lists.parse: Reading linelist ---
   /Users/ncrighton/Code/Repo/linetools/build/lib.macosx-10.5-x86_64-3.4/linetools/data/lines/morton03_table2.fits.gz
linetools.lists.parse: Reading linelist ---
   /Users/ncrighton/Code/Repo/linetools/build/lib.macosx-10.5-x86_64-3.4/linetools/data/lines/morton00_table2.fits.gz
linetools.lists.parse: Reading linelist ---
   /Users/ncrighton/Code/Repo/linetools/build/lib.macosx-10.5-x86_64-3.4/linetools/data/lines/verner96_tab1.fits.gz
linetools.lists.parse: Reading linelist ---
   /Users/ncrighton/Code/Repo/linetools/build/lib.macosx-10.5-x86_64-3.4/linetools/data/lines/verner94_tab6.fits
linetools.lists.parse: Reading linelist ---
   /Users/ncrighton/Code/Repo/linetools/build/lib.macosx-10.5-x86_64-3.4/linetools/data/lines/EUV_lines.ascii
read_sets: Using set file --
  /Users/ncrighton/Code/Repo/linetools/build/lib.macosx-10.5-x86_64-3.4/linetools/lists/sets/llist_v1.0.ascii

# HILyman system
HIsys = lt_absys.LymanAbsSystem.from_components([abscomp])
print(HIsys)
print(HIsys._components)

<LymanAbsSystem: name=J081227.432-122555.56_z2.929 type=HILyman, 08:12:27.432 -12:25:55.56, z=2.92939, NHI=0>
[<AbsComponent: 08:12:27.432 -12:25:55.56, Name=HI_z2.92939, Zion=(1,1), Ej=0 1 / cm, z=2.92939, vlim=-300 km / s,300 km / s>]
Multiple components 
# SiII
SiIItrans = ['SiII 1260', 'SiII 1304', 'SiII 1526', 'SiII 1808']
abslines = []
for trans in SiIItrans:
    iline = AbsLine(trans)
    iline.attrib['z'] = 2.92939
    iline.analy['vlim'] = [-250.,80.]*u.km/u.s
    abslines.append(iline)
#
SiII_comp = AbsComponent.from_abslines(abslines)
SiII_comp.coord = radec

# Generic
imp.reload(lt_absys)
LLSsys = lt_absys.GenericAbsSystem.from_components([abscomp,SiII_comp])
print(LLSsys)
print(LLSsys._components)

<GenericAbsSystem: name=Foo type=Generic, 08:12:27.432 -12:25:55.56, z=2.92939, NHI=0>
[<AbsComponent: 08:12:27.432 -12:25:55.56, Name=HI_z2.92939, Zion=(1,1), Ej=0 1 / cm, z=2.92939, vlim=-300 km / s,300 km / s>, <AbsComponent: 08:12:27.432 -12:25:55.56, Name=SiII_z2.92939, Zion=(14,2), Ej=0 1 / cm, z=2.92939, vlim=-250 km / s,80 km / s>]
Methods
List of AbsLines 
lines = LLSsys.list_of_abslines()
lines

[<AbsLine: HI 1215, wrest=1215.6700 Angstrom>,
 <AbsLine: HI 1025, wrest=1025.7222 Angstrom>,
 <AbsLine: SiII 1260, wrest=1260.4221 Angstrom>,
 <AbsLine: SiII 1304, wrest=1304.3702 Angstrom>,
 <AbsLine: SiII 1526, wrest=1526.7070 Angstrom>,
 <AbsLine: SiII 1808, wrest=1808.0129 Angstrom>]
Single Line 
lyb = LLSsys.get_absline('HI 1025')
lyb

<AbsLine: HI 1025, wrest=1025.7222 Angstrom>

lyb = LLSsys.get_absline(1025.72*u.AA)
lyb

<AbsLine: HI 1025, wrest=1025.7222 Angstrom>

lyb.wrest
\[1025.7222 \; \mathrm{\mathring{A}}\]

Overview

This Class is designed to organize and analyze an absorption system. This is generally constructed of one or more AbsComponent Class. The base class is abstract, i.e. one must instantiate one of its flavors (e.g. HILyman, MgII, LLS, DLA).

By definition, an AbsSystem is a unique collection of absorption components. It is specified by:

Property Variable Type Description
RA, DEC radec tuple or coord RA,DEC in deg or astropy.coordinate
Redshift z float absorption redshift
Velocity limits vlim Quantity array -/+ velocity limits of the system

Instantiation

The AbsSystem Class may be instantiated in a few ways. The default sets the properties listed above:

gensys = GenericAbsSystem((15.23*u.deg,-23.22*u.deg), 1.244, [-500,500]*u.km/u.s, NHI=16.)

More commonly, one will instantiate with one or more AbsComponent objects:

# HI Lya, Lyb
radec = SkyCoord(ra=123.1143*u.deg, dec=-12.4321*u.deg)
lya = AbsLine(1215.670*u.AA, z=2.92939)
lya.limits.set([-300.,300.]*u.km/u.s)  # vlim
lyb = AbsLine(1025.7222*u.AA, z=lya.attrib['z'])
lyb.limits.set([-300.,300.]*u.km/u.s)  # vlim
abscomp = AbsComponent.from_abslines([lya,lyb])
abscomp.coord = radec
# Finish
HIsys = LymanAbsSystem.from_components([abscomp])

One may also instantiate from a dict, usually read from the hard-drive:

abscomp = AbsSystem.from_dict(idict)

Attributes

Sub Classes

Generic

A catch-all subclass for AbsSystem. More options are provided in pyigm.

Plots

Methods

AbsLines

There are a few methods related to the AbsLine objects within an AbsSystem. One can generate a list of all the AbsLine objects with:

lines = abssys.list_of_abslines()

One can retrieve one or more AbsLine objects matching the name or rest-wavelength of a transition, e.g.

lyb = abssys.get_absline('HI 1025')
# or
lyb = abssys.get_absline(1025.72*u.AA)  # Nearest 0.01 A is required
Components

get_component grabs components matching an input where the input is either a tuple of (Z, ion) or an AbsLine:

SiII = gensys.get_component((14,2))

update_component_colm synthesizes and updates the column densities for the components.:

gensys.update_component_colm()
ionN

Fill the _ionN attribute with a QTable of column densities. These are derived from the components

abssys.fill_ionN()
print(abssys._ionN)

Output

One may generate a dict of the key properties of the AbsSystem with the to_dict() method:

odict = HIsys.to_dict()

This dict is required to be JSON compatible.

RelAbund Class

Notebooks

Overview

This class packages the relative abundances of an object, typically an AbsSystem.

Instantation

Init

One can instantiate via the init and then fill the data dict. This is a bit cumbersome and not especially recommended. But here is an example.

To begin, make a new class instance:

XY = RelAbund()
Loading abundances from Asplund2009

Then load data into the data dict. Here is an example:

XY._data = {6: dict(flag=1, XH=-1., sigXH=0.2, sig=0.05),
         8: dict(flag=2, XH=-1.4, sigXH=0.25, sig=0.05),
         14: dict(flag=1, XH=-1.1, sigXH=0.25, sig=0.05),
         26: dict(flag=1, XH=-1.4, sigXH=0.25, sig=0.05),
         32: dict(flag=3, XH=-0.8, sigXH=0.25, sig=0.05),
         }

The flag value indicate the type of measurement:

Flag Description
1 Standard value (and error)
2 Lower limit (e.g. saturated line)
3 Upper limit (e.g. blend or non-detection)
Ionic Column Table

More frequent usage will be to instantiate using an input table of column density measurements, e.g.:

dla.XY = RelAbund.from_ionclm_table((1,dla.NHI, dla.sig_NHI[0]), dla._ionN)

See pyigm DLA abund Notebook for more.

By Hand

For quick and dirty abundance calculations, you may find the from_pair method useful:

Usage

You may grab the data for any element with item syntax:

CH = XY[6]
{'flag': 1, 'sig': 0.2, 'val': -1.0}
CH = XY['C']

Element ratios can be accessed by providing a tuple of atomic number or element name:

SiFe = XY['Si', 'Fe']
{'flag': 1, 'sig': 0.070710678118654766, 'val': 0.2999999999999998}

You may generate an astropy Table of the X/Y values:

tbl = XY.table()  # For X/H
tbl = XY.table('Fe')  # For X/Fe

SolarAbund Class

Notebooks

Examples with the SolarAbund Class (v1.1)

Download this notebook.

# import
from linetools.abund import solar as labsol
Init
sol = labsol.SolarAbund()

Loading abundances from Asplund2009
Abundances are relative by number on a logarithmic scale with H=12

sol

<SolarAbund: Asplund2009>
Usage
# Simple calls
print(sol['C'])
print(sol[6])

8.43
8.43

# Ratio
print(sol.get_ratio('C/Fe'))

0.98

# Multiple calls
print(sol[6,12,'S'])

[ 8.43  7.53  7.15]
Bits and pieces
from linetools.abund import ions as laions

# Ion name
laions.ion_name((6,2))

'CII'

# Name to ion
laions.name_ion('CII')

(6, 2)

from linetools.abund.elements import ELEMENTS

ele = ELEMENTS['C']

ele.eleconfig_dict

{(1, 's'): 2, (2, 'p'): 2, (2, 's'): 2}

Overview

This class provides access to element abundance ratios measured in (or nearby) the Sun.

To access the abundances, make a new class instance:

>>> from linetools.abund import solar as labsol
>>> sol = labsol.SolarAbund()
Loading abundances from Asplund2009
Abundances are relative by number on a logarithmic scale with H=12

Then select the element you want by either its name or atomic number:

>>> sol['C']
8.4299999999999997
>>> sol[6]
8.4299999999999997

Currently the abundances from Asplund et al. 2009 are available, and in future more references will be included.

Multiple elements can also be selected:

>>> sol['C', 'O']
array([ 8.43,  8.69])

Element ratios can be accessed using the get_ratio method:

>>> sol.get_ratio('C/Fe')
0.97999999999999954

LineList Class

Overview

This class organizes information about atomic and molecular transition lines (e.g. HI Lya, CIV 1548, Hydrogen Balmer series) observed in astrophysical environments.

The following lists are currently avaliable:

  • ‘ISM’ : “All” ISM lines (can be overwhelming!)
  • ‘Strong’ : Strong ISM lines (most common absorption line transitions observed)
  • ‘HI’ : HI Lyman series
  • ‘H2’ : H2 (Lyman-Werner)
  • ‘CO’ : CO UV band-heads
  • ‘EUV’ : Extreme UV lines
  • ‘Galaxy’ : Lines typically identified in galaxy spectra

Instantiation

The LineList Class may be instantiated using one of the keys in the list above:

from linetools.lists.linelist import LineList
hi = LineList('HI')
#
euv = LineList('EUV')

hi, for example, contains only HI Lyman series transitions (e.g. HI Lya), and euv contains both HI Lyman series and extreme UV metal transitions (e.g. HI Lyb, NeVIII, MgX).

Accessing single transitions

We can now easily access atomic information regarding individual transitions either by the rest-frame wavelength:

wrest = 1215.67 * u.AA  # HI Lya
hi[wrest]

or by the name convention within linetools, which in the case of HI Lya is HI 1215:

name = 'HI 1215'   #  We adopt the convention of *not* rounding in the name
hi[name]

Both cases will provide the following dictionary:

{'A': <Quantity 626500000.0 1 / s>,   # Einstein coefficient
'Am': 0,                              # Mass number (often written as "A"; only used for D)
'Ej': <Quantity 0.0 1 / cm>,          # Energy of lower level (relative to ground state)
'Ek': <Quantity 2259.163 1 / cm>,     # Energy of upper level (relative to ground state)
'Ex': <Quantity 0.0 1 / cm>,          # Excitation energy (cm^-1)
'Jj': 0.0,                            # Tot ang mom (z projection) of lower state (or rotation level)
'Jk': 0.0,                            # Tot ang mom (z projection) of upper state (or rotation level)
'Ref': 'Morton2003',                  # References
'Z': 1,                               # Atomic number (for atoms)
'col0': masked,                       # (Reserved)
'col6': masked,                       # (Reserved)
'el': 0,                              # Electronic transition (2=Lyman (B-X), 3=Werner (C-X))
'f': 0.41639999999999999,             # Oscillator strength
'gamma': <Quantity 626500000.0 1 / s>,# Sum of A
'gj': 2,                              # Lower statistical weight (2J+1)
'gk': 6,                              # Upper statistical weight (2J+1)
'group': 1,                           # Flag for grouping
'ion': 1,                             # Ionic state (1=Neutral)
'mol': '',                            # Molecular name (H2, HD, CO, C13O)
'name': 'HI 1215',                    # Name
'nj': 0,                              # Orbital level of lower state (or vibrational level)
'nk': 0,                              # Orbital level of upper state (or vibrational level)
'wrest': <Quantity 1215.67 Angstrom>} # Rest Wavelength (Quantity)

which summarizes the most important atomic information of HI Lya transition, including the reference where these values come from (i.e., Morton2003). One can therefore access any of these by calling its dictionary keywords:

hi['HI 1215']['wrest']
<Quantity 1215.67 Angstrom>

is the rest-frame wavelength of the HI Lya transition. Similarly,:

euv['NeVIII 780']['f']
0.050500001758337021

is the oscillator strength of the NeVIII 780 transition.

Methods

subset_lines()

This method provides a way to define a subset of lines drawn from the original` LineList` object. Consider that for some reason you may want only HI Lya and Lyb in your LineList, then you can achieve this by:

hi = LineList('HI')
hi = hi.subset_lines(['HI 1215', 'HI 1025'])

Which has only those two transitions loaded.

You may also want to use rest-frame wavelength to define a subset, for instance:

ism = LineList('ISM')
lines = [2796.3543, 2803.5315, 1548.195, 1550.77] * u.AA
ism = ism.subset_lines(lines)
print(ism)
<LineList: ISM; 4 transitions>

selects only those four transitions of MgII and CIV. In order to avoid loading the LineList('ISM') again, you can use the keyword reset_data in subset_lines() to make another arbitrarily different subset of lines from the original LineList:

lines = ['HI 1215', 'HI 1025']
ism = ism.subset_lines(lines, reset_data=True)
print(ism)
<LineList: ISM; 2 transitions>

which now has only HI Lya and Lyb.

Finally, if you want the transitions to be sorted by rest-frame wavelength you can use the optional keyword sort:

lines = [2796.3543, 2803.5315, 1548.195, 1550.77] * u.AA
ism = ism.subset_lines(lines, reset_data=True, sort=True)
ism._data['wrest']
<Quantity [ 1548.195 , 1550.77  , 2796.3543, 2803.5315] Angstrom>
set_lines()

Another way to reset the LineList to its original form is by using set_lines(). Following the previous example, we have a ism Linelist with only 4 transitions:

print(ism._data['name'])
   name
---------
CIV 1548
CIV 1550
MgII 2796
MgII 2803

print(ism)
<LineList: ISM; 4 transitions>

ism.set_lines()
print(ism)
<LineList: ISM; 412 transitions>

Give us the original ism LineList with 412 unique transitions.

You may also want to use rest-frame wavelength to define a subset, for instance:

ism = LineList('ISM')
sub_lines = [2796.3543, 2803.5315, 1548.195, 1550.77] * u.AA
civ_mgii = ism.subset(sub_lines)
all_transitions()

Sometimes it may be useful to know all the transitions associated to a given ion species. This can be achieved by the all_transitions() method:

ism = LineList('ISM')
mgii = ism.all_transitions('MgII')

Which give us the information of all the 6 transitions of MgII:

print(mgii)
     A       el  nj  nk group    name       Ek    ...  Jk  Z   gk  gj    gamma    col0 col6
    1 / s                                  1 / cm  ...                    1 / s
----------- --- --- --- ----- --------- --------- ... --- --- --- --- ----------- ---- ----
  2350000.0   0   0   0     1 MgII 1025  97468.92 ... 0.0  12   4   2   2350000.0   --   --
  2480000.0   0   0   0     1 MgII 1026  97455.12 ... 0.0  12   2   2   2480000.0   --   --
  1370000.0   0   0   0     1 MgII 1239  80650.02 ... 0.0  12   4   2   1370000.0   --   --
  1540000.0   0   0   0     1 MgII 1240   80619.5 ... 0.0  12   2   2   1540000.0   --   --
262500000.0   0   0   0     1 MgII 2796 35760.848 ... 0.0  12   4   2 262500000.0   --   --
259500000.0   0   0   0     1 MgII 2803 35669.298 ... 0.0  12   2   2 259500000.0   --   --

In this case mgii is a QTable because more than 1 transition was found. In cases were only 1 transition exists, the output of all_transitions() is a dictionary with the same keywords as the columns of ism._data QTable:

ciii = ism.all_transitions('CIII')
type(ciii)
dict
print(ciii)
{'A': <Quantity 1760000000.0 1 / s>,
'Am': 0,
'Ej': <Quantity 0.0 1 / cm>,
'Ek': <Quantity 2352.04 1 / cm>,
'Ex': <Quantity 0.0 1 / cm>,
'Jj': 0.0,
'Jk': 0.0,
'Ref': 'Morton2003',
'Z': 6,
'col0': masked,
'col6': masked,
'el': 0,
'f': 0.75700000000000001,
'gamma': <Quantity 1760000000.0 1 / s>,
'gj': 1,
'gk': 3,
'group': 1,
'ion': 3,
'mol': '',
'name': 'CIII 977',
'nj': 0,
'nk': 0,
'wrest': <Quantity 977.0201 Angstrom>}

You can also use a rest-frame wavelength to identify the ion species of interest:

wrest = 1260.4221 * u.AA
si2 = ism.all_transitions(wrest)
print(si2['name', 'wrest', 'f'])
   name     wrest          f
           Angstrom
--------- --------- ---------------
SiII 889  889.7228 0.0434000007808
SiII 989  989.8731           0.171
SiII 1020 1020.6989          0.0168
SiII 1190 1190.4158           0.292
SiII 1193 1193.2897           0.582
SiII 1260 1260.4221            1.18
SiII 1304 1304.3702          0.0863
SiII 1526  1526.707           0.127
SiII 1808 1808.0129         0.00208
SiII 2335  2335.123        4.25e-06

For the purposes of all_transitions, it does not matter which transition of a given ion species you choose, it will still retrieve the same answer, e.g.:

hi = ism.all_transitions('HI 1215')
hi = ism.all_transitions('HI 1025')
hi = ism.all_transitions(972.5367 * u.AA)
hi = ism.all_transitions('HI')

are all equivalent. Note that in the last example we only used the root name of the transition (i.e. the string before the blank space, 'HI'), so no prior knowledge of the linetools naming convention is needed.

strongest_transitions()

Sometimes it is useful to know the strongest transition for an ion in the LineList within some wavelength range. strongest_transitions() gives the strongest n_max transitions of a given ion between a wavelength range, sorted by relative strength (defined as the product of its rest-frame wavelength wrest and oscillator strength f):

wvlims = [1000, 3000] * u.AA
line = 'SiII'
si2_strong = ism.strongest_transitions(line, wvlims, n_max=4)
print(si2_strong['name'])
   name
---------
SiII 1260
SiII 1193
SiII 1190
SiII 1526

The syntax is the same as for all_transitions(). Note that you will get the same result if you use line='SiII', line='SiII 1190', line='SiII 889', or line=889.7228*u.AA. By default n_max=3. Depending on the wavelength range, however, the output may vary:

wvlims = [500, 1100] * u.AA
line = 'SiII 1260'
si2_strong = ism.strongest_transitions(line, wvlims, n_max=4)
print(si2_strong['name'])
   name
---------
SiII 989
SiII 889
SiII 1020

Note that despite n_max=4 we have only retrieved the 3 transitions satisfying the criteria of belonging to wvlims = [500, 1100] * u.AA. Again, note that even though SiII 1260 is out of wvlims range, it can still be used to identify that you are interested in the SiII ion species.

If you would like to retrieve all the transitions in a given wvlims regardless of its relative strength, you can set n_max=None.

Following the convention within LineList, if only 1 transition is retrieved, the output of strongest_transitions() is a dictionary; if more than 1 transition are retrieved the output is a QTable. If no transition exist the output is None.

available_transitions()

Sometimes it may be useful to know what are the available transition in a given wavelength range found in the LineList regardless of the ion species. This is particularly the case when someone is trying to identify unknown emission/absorption lines in a spectrum. Let us then illustrate the use of this method with an example. Imagine that you have an observed spectrum covering the following wavelength range:

wvlims = [3500,5000] * u.AA

Let us now imagine that we are interested in a particular redshift, say z=0.67. Then, we can do:

z = 0.67
transitions = ism.available_transitions(wvlims/(1+z), n_max_tuple=None, min_strength=0.)
print(len(transitions))
33

Will give the 33 transitions available that could correspond to having z=0.67 in the form of a QTable. The output is sorted by strength of the strongest available transition per ion species, and strength is defined as log10(wrest * fosc * abundance), where abundance is that of the solar composition given by Asplund2009. As optional keyword parameters one can specify a minimum strength as min_strength, so transitions below this value are omitted, e.g.:

transitions = ism.available_transitions(wvlims/(1+z), n_max_tuple=None, min_strength=10.5)
print(len(transitions))
3

Which correspond to MgI 2852, MgII 2796 and MgII 2803. Note than this method does not correct for ionization state. Similarly, one can also specify the maximum number of transitions per ion species tuple using the optional keyword parameter n_max_tuple, e.g.:

transitions = ism.available_transitions(wvlims/(1+z), n_max_tuple=1, min_strength=0.)
print(transitions['name'])
    name
-----------
MgI 2852
MgII 2796
FeII 2382
FeII* 2396b
MnII 2576
VII  2683
...

Which for the case of MgII only retrieves 'MgII 2796'. Again, following the convention within LineList, if only 1 transition is retrieved, the output of available_transitions() is a dictionary; if more than 1 transition are retrieved the output is a QTable. If no transition exist satisfying the criteria the output is None.

SpectralLine Class

Notebooks

Overview

This Class is designed to organize and analyze an spectral line, either emission or absorption. This is an abstract base class, i.e. one must instantiate one of its subclasses (e.g. AbsLine).

Sub Classes

The primary children of SpectralLine are AbsLine Class and EmissionLine (to be implemented). See their documentation for a description of instantiation and additional attributes.

Attributes

The base attributes for the SpectralLine class are:

Property Variable Type Description
RA, Dec attrib[‘coord’] Coord astropy.coordinate
Redshift attrib[‘z’] float Reference redshift
Redshift sigma attrib[‘sig_z’] float Reference redshift uncertainty
Velocity attrib[‘v’] Quantity line velocity relative to its redshift
Velocity sigma attrib[‘sig_v’] Quantity 1 sigma uncertainty in the velocity
Equivalent Width attrib[‘EW’] Quantity Equivalent width
EW sigma attrib[‘sig_EW’] Quantity 1 sigma uncertainty in EW
EW flag attrib[‘flag_EW’] int Equivalent width flag
Limits limits LineLimits The limits of the line in redshift, velocity (w/r to its redshift) and observed wavelength.

Analysis

It is common that one wishes to associate a line with a spectrum to perform a range of analyses. This is accomplished through:

spline.analy['spec'] = sp

where sp is an XSpectrum1D Class object.

Methods

cut_spec

Provide a spectrum has been associated to the line (see Analysis): then this method returns the portion of the spectrum surrounding the line. The limits are specified in the LineLimits class held in the attribute limits, usually either with observed wavelengths or velocities relative to the line’s redshift. The code returns the flux, error array, and a dict containing the wavelength and velocity arrays.

spline.limits.set([-300., 300.]*u.km/u.s) # vlim
fx, sig, wv_dict = spline.cut_spec()
ismatch

Check whether two Lines are equal rest wavelength (to 1e-6 AA), whether they have common RA/DEC to 0.1” (if provided), and whether they are from the same ion:

print specline.ismatch(another_line)
measure_ew

Measure the observed Equivalent width for a SpectralLine. Following absorption-line convention, absorption will have a positive value and emission will have a negative value.

To perform the calculation, the line must be associated to a spectrum (see Analysis_) and either wvlim or vlim must be specified. When executed, the EW and sig_EW attibutes are filled:

specline.measure_ew()
measure_kin

Measure kinematic characteristics of an AbsLine. To perform the calculation, the line must be associated to a spectrum (see Analysis_) and vlim must be specified. When executed, the ‘kin’ attribute is filled with a dict of measurements. Default set of measurements are the v90, fedg, and fmm statistics of Prochaska & Wolfe 1997:

specline.measure_kin()
measure_restew

Measure the rest-frame Equivalent width of a SpectralLine. See measure_ew for details.

to_dict

Convert the Class to a JSON-friendly dict that might be easily written to the disk, e.g.:

adict = specline.to_dict()
with io.open(outfil, 'w', encoding='utf-8') as f:
   f.write(unicode(json.dumps(tdict, sort_keys=True,
      indent=4, separators=(',', ': '))))

Utilities

There are several utilites related to spectral lines. These are located in the line_utils module.

parse_speclines

Given a list of SpectralLines and desired property (key), this method returns a list or array of the values.:

from linetools import line_utils
array_of_values = line_utils.parse_speclines(list_of_speclines, mk_array=True)
transtable_from_speclines

Given a list of SpectralLines, this method returns a Table of a subset of the properties (e.g. wavelength, name, EW).:

trans_tbl = line_utils.transtable_from_speclines(list_of_speclines)

XSpectrum1D Class

Overview

XSpectrum1D describes a 1-d spectrum, which usually consists of a wavelength, flux and flux uncertainty array. For absorption-line analysis, it also often contains a continuum array. The data are held in a masked numpy array which may contain multiple spectra. By default pixels on the edges of the spectrum with an input error array having values <=0 are masked. It is important to appreciate this masking. It does mean that you will not view, print, analyze, etc. pixels that have been masked.

Attributes

The main attributes of XSpectrum1D are wavelength, flux and sig. Let’s begin by creating a spectrum using the from_tuple method:

>>> from linetools.spectra.xspectrum1d import XSpectrum1D
>>> import numpy as np
>>> wa = np.arange(3000, 7000.1, 0.5)
>>> fl = np.ones_like(wa)
>>> sig = np.ones_like(fl) * 0.1
>>> sp = XSpectrum1D.from_tuple((wa, fl, sig), verbose=False)
>>> sp.wavelength # doctest: +SKIP
<Quantity [ 3000. , 3000.5, 3001. ,..., 6999. , 6999.5, 7000. ] Angstrom>
>>> sp.flux
<Quantity [ 1., 1., 1.,...,  1., 1., 1.]>
>>> sp.sig # doctest: +SKIP
<Quantity [ 1., 1., 1.,...,  1., 1., 1.]>

Note that all three arrays have units. If you don’t specify a unit when you create an new XSpectrum1D instance, Angstroms are assumed for wavelength and dimensionless_unscaled for flux. The 1-sigma uncertainty is always assumed to have the same units as the flux. All of these are specified in the sp.units dict.

If one loads multiple 1D spectra (e.g. a brick of data from DESI or a set of spectrum from igmspec), the selected spectrum is given by the spec.select index.

All of the values are stored in the masked spec.data numpy array with columns wave, flux, sig, and co (the latter is for a continuum).

Init

Reading

Read spectra from a file using XSpectrum1D.from_file, which uses the same syntax as readspec. See below for a complete listing of permitted file formats.

The easiest way to create a new spectrum from a set of data arrays for a single spectrum is to use sp.from_tuple as shown above. Here are a series of example calls to generate the class:

sp = XSpectrum1D.from_file('PH957_f.fits')      # From a FITS file
sp = XSpectrum1D.from_file('q0002m422.txt.gz')  # From an ASCII table
sp = xspec1.copy()                              # From an XSpectrum1D object
sp = XSpectrum1D.from_tuple((wa, fl, sig), verbose=False)
Masking

The guts of XSpectrum1D is a ndarray array named data which contains the wave, flux, sig, etc. values. This is a masked array which is convenient for many applications. If you wish to view/analyze all pixels in your spectrum including those with 0 or NAN sig values, then disable the mask when creating the object or by using the unnmask() method:

sp = XSpectrum1D.from_tuple((wa, fl, sig), masking='none')
sp = XSpectrum1D.from_file('PH957_f.fits')
sp.unmask()

Methods

Writing

There are a number of methods to write a file, e.g. sp.write_to_fits. FITS files are preferable because they are generally faster to read and write, require less space, and are generally easier for other software to read. Another option is an HDF5 file which better preserves the data format of XSpectrum1D. Here are some examples:

sp.write_to_fits('QSO.fits')            # Standard FITS file
sp.write('QSO.fits')                    # Same
sp.write('QSO.fits', FITS_TABLE=True)   # Binary FITS table
sp.write_to_hdf5('QSO.hdf5')            # HDF5 file
sp.write('QSO.hdf5')                    # Same
sp.write_to_ascii('QSO.ascii')          # ASCII (heaven forbid)
sp.write('QSO.ascii')                   # Same

One can collate a list of XSpectrum1D objects into one with collate:

sp1 = XSpectrum1D.from_file('PH957_f.fits')
sp2 = XSpectrum1D.from_file('q0002m422.txt.gz')
sp = linetools.spectra.utils.collate([sp1,sp2])
Plotting

sp.plot() plots the spectrum, which you can then navigate around using the same keys as lt_xspec (as well as the usual matplotlib navigation tools). Note: if you are using MacOSX then you will probably need to change your backend from macosx to TkAgg in the matplotlibrc file.

Rebinning

rebin rebins the spectrum to an arbitrary input wavelength array. Flux is conserved. If do_sig=True, the error array is rebinned as well and a crude attempt is made to conserve S/N. Generally, neighboring pixels will be correlated:

newspec = sp.rebin(new_wv, do_sig=True)

If the XSpectrum1D object containts multiple spectra, you can rebin all of them to the new wavelength array as well:

newspec = sp.rebin(new_wv, do_sig=True, all=True)
Continuum fitting

fit_continuum enables you to interactively fit a continuum to the spectrum. Currently it’s optimised to estimate the continuum for high-resolution quasar spectra, but it should be applicable to any spectrum with a slowly varying continuum level and narrow absorption features. Once a continuum has been fitted, it can be accessed using the co attribute. The spectrum can also be normalised (i.e the flux values returned by spec.flux are divided by the continuum) with the normalize method. This also sets spec.normed to True.

Finally, you can apply small variations to the continuum anchor points with perturb_continuum to see how changes in the continuum level affect your analysis.

Smoothing

There are several algorithms included that smooth the input spectrum and return a new XSpectrum1D. These are box_smooth, gauss_smooth, and ivar_smooth.

Other methods

You can join one XSpectrum1D instance with another overlapping spectrum using splice. pix_minmax finds the pixel indices corresponding to a wavelength or velocity range, and add_noise adds noise to the spectrum. We have also implemented a method that estimates a local average signal-to-noise ratio at a given observed wavelength (get_local_s2n), which is capable of masking out pixels that are below a flux threshold (useful for excluding strong absorption features from the calculation). For a complete list of all the available methods, see the API: XSpectrum1D.

Multi-spec methods

See Multi XSpectrum1D for more.

File Formats Read

Below is a table of the types of spectra files that can be read by readspec. If your file cannot be read, please open an issue on the linetools issue tracker.

Description Instruments
simple 1D FITS files ESI, HIRES, etc.
binary FITS table from LowRedux LRIS,Kast,etc.
multi-extension 1D FITS files from LowRedux LRIS,Kast,etc.
binary FITS tables from many other sources COS, SDSS, etc.
multi-extension binary FITS tables from PYPIT LRIS,Kast,etc.
brick files (2D images: flux, ivar; 1D image: wavelength) DESI
UVES_popler output files UVES

Graphical User Interfaces (GUIs)

GUIs

Notebooks

Interactive continuum fitting

Download this notebook.

# suppress warnings for these examples
import warnings
warnings.filterwarnings('ignore')

import imp
prefix = imp.find_module('linetools')[1] + '/spectra/tests/files/'
import linetools.spectra.xspectrum1d as lsx
spec = lsx.XSpectrum1D.from_file(prefix + 'q0002m422.txt.gz')
# keep the old continuum to compare later on
co_old = spec.co.copy()

%pylab

Using matplotlib backend: TkAgg
Populating the interactive namespace from numpy and matplotlib

# now fit the continuum interactively. We say we're fitting a QSO,
# so it can make intelligent guesses about where to put the spline
# points that define the continuum.
spec.fit_continuum(kind='QSO', redshift=2.76)

# now you can interactively tweak these spline points, adding or
# removing them as necessary. Once you're finished, press 'q' to
# close the window.

knots file exists, use this? (y) n

i,o          Zoom in/out x limits
y            Zoom out y limits
Y            Guess y limits
t,b          Set y top/bottom limit
l,r          Set left/right x limit
[,]          Pan left/right
w            Plot the whole spectrum

S,U          Smooth/unsmooth spectrum


i,o          Zoom in/out x limits
y            Zoom out y limits
Y            Guess y limits
t,b          Set y top/bottom limit
l,r          Set left/right x limit
[,]          Pan left/right
w            Plot the whole spectrum

S,U          Smooth/unsmooth spectrum

a        : add a new spline knot
A        : add a new spline knot, and use a flux median to guess y position
+        : double the number of spline knots
_        : halve the number of spline knots
d        : delete the nearest knot
m        : move the nearest knot
M        : move the nearest knot, and use a flux median to guess y position

q        : quit

Updating continuum

# the New continuum is now saved in spec.co, and the spline knots are in
# spec.meta['contpoints']
#
# Let's compare the old and new continuum
plt.figure()
wa = spec.dispersion.value
plt.plot(wa, co_old)
plt.plot(wa, spec.co)

[<matplotlib.lines.Line2D at 0x10c8df978>]

co_old2 = spec.co.copy()

# we can also tweak a small section of the continuum without affecting the whole spectrum.
spec.fit_continuum(wlim=(5000, 5100))

knots file exists, use this? (y) n

i,o          Zoom in/out x limits
y            Zoom out y limits
Y            Guess y limits
t,b          Set y top/bottom limit
l,r          Set left/right x limit
[,]          Pan left/right
w            Plot the whole spectrum

S,U          Smooth/unsmooth spectrum


i,o          Zoom in/out x limits
y            Zoom out y limits
Y            Guess y limits
t,b          Set y top/bottom limit
l,r          Set left/right x limit
[,]          Pan left/right
w            Plot the whole spectrum

S,U          Smooth/unsmooth spectrum

a        : add a new spline knot
A        : add a new spline knot, and use a flux median to guess y position
+        : double the number of spline knots
_        : halve the number of spline knots
d        : delete the nearest knot
m        : move the nearest knot
M        : move the nearest knot, and use a flux median to guess y position

q        : quit

Updating continuum

# check it works without a predefined continuum
spec = lsx.XSpectrum1D.from_file(prefix + 'q0002m422.txt.gz')
spec.co = None
spec.fit_continuum(kind='QSO', redshift=2.76)

knots file exists, use this? (y) n

i,o          Zoom in/out x limits
y            Zoom out y limits
Y            Guess y limits
t,b          Set y top/bottom limit
l,r          Set left/right x limit
[,]          Pan left/right
w            Plot the whole spectrum

S,U          Smooth/unsmooth spectrum


i,o          Zoom in/out x limits
y            Zoom out y limits
Y            Guess y limits
t,b          Set y top/bottom limit
l,r          Set left/right x limit
[,]          Pan left/right
w            Plot the whole spectrum

S,U          Smooth/unsmooth spectrum

a        : add a new spline knot
A        : add a new spline knot, and use a flux median to guess y position
+        : double the number of spline knots
_        : halve the number of spline knots
d        : delete the nearest knot
m        : move the nearest knot
M        : move the nearest knot, and use a flux median to guess y position

q        : quit

Updating continuum
xspec Documentation

Download this notebook.

Download this notebook.

This ipython Notebook is intended to provide documentation for the linetools GUI named XSpecGUI.

Enjoy and feel free to suggest edits/additions, etc.

Here is a screenshot of the XSpecGUI in action:

from IPython.display import Image
Image(filename="images/xspec_example.png")
_images/xspecgui_1_0.png

The example spectrum file used below is part of the linetools package.

import imp
lt_path = imp.find_module('linetools')[1]
spec_fil = lt_path+'/spectra/tests/files/PH957_f.fits'
Before Launching the GUI

If you are a Mac user, we highly recommend that you set your matplotlib backend from MacOSX to TkAgg (or another option, see backends).

Launching the GUI
From within ipython or equivalent 
from linetools.guis import xspecgui as ltxsg

import imp; imp.reload(ltxsg)
ltxsg.main(spec_fil)
Overlaying Line Lists

You can overlay a series of vertical lines at standard spectral lines at any given redshift.

Setting the Line List

You must choose a line-list by clicking one.

Setting the redshift
  • Type one in
  • RMB on a spectral feature (Ctrl-click on Emulated 3-button on Macs)
    • Choose the rest wavelength
Marking Doublets
  • C – CIV
  • M – MgII
  • X – OVI
  • 4 – SiIV
  • 8 – NeVIII
  • B – Lyb/Lya
Velocity plot (Coming Soon)

Once a line list and redshift are set, type ‘v’ to launch a Velocity Plot GUI.

Simple Analysis
Gaussian Fit

You can fit a Gaussian to any single feature in the spectrum as follows: 1. Click “G” at the continuum at one edge of the feature 1. And then another “G” at the other edge (also at the continuum) 1. A simple Gaussian is fit and reported.

Equivalent Width

You can measure the rest EW of a spectral feature as follows: 1. Click “E” at the continuum at one edge of the feature 1. And then another “E” at the other edge (also at the continuum) 1. A simple boxcar integration is performed and reported.

Apparent Column Density

You can measure the apparent column via AODM as follows: 1. Click “N” at the continuum at one edge of the feature 1. And then another “EN” at the other edge (also at the continuum) 1. A simple AODM integration is performed and reported.

Ly\(\alpha\) Lines
  • “D” - Plot a DLA with \(N_{\rm HI} = 10^{20.3} \rm cm^{-2}\)
  • “R” - Plot a SLLS with \(N_{\rm HI} = 10^{19} \rm cm^{-2}\)
xabssys Documentation

This ipython Notebook is intended to provide documentation for the linetools GUI named XAbsSysGui.

Enjoy and feel free to suggest edits/additions, etc.

Before Launching the GUI

If you are a Mac user, we highly recommend that you set your matplotlib backend from MacOSX to TkAgg (or another option, see backends).

Launching the GUI
From within ipython or equivalent

Not yet implemented

Modifying Absorption lines
Limits and Blends

Overview

There are several GUIs included with linetools, primarily for simple spectral inspection and analysis.

We caution that it is difficult (essentially impossible) to generate unit tests for these GUIs. As such, they are far from bug free and may crash unexpectedly. Buyer beware!

Continuum fitting

This enables interactive fitting of the unabsorbed continuum for a spectrum. A series of ‘knot’ positions are estimated, and these are then joined with a spline to produce a continuum. Spline points can be interactively added, deleted or moved to improve the continuum. See the notebook for an example.

XSpecGUI

This enables visual inspection of a spectrum. Simple analysis (e.g. equivalent width measurements) may also be performed. See the notebook for details.

XAbsSysGui

This shows a velocity (stack) plot of the absorption lines from an input absorption line system. The user can then modify the velocity limits that would be used for subsequent analysis, flag bad lines, blends, set limits, etc. The modified absoprtion system is then written to the hard-drive as a JSON file.

Command line tools

linetool Scripts

There are a number of scripts, many of which are GUIs, provided with linetools. As regards the GUIs we warn again that Mac users will need to set their matplotlib to something other than MacOSX. See backends.

lt_xspec

Launch a QT Gui to examine an input spectrum.

For example:

lt_xspec filename.fits
lt_xspec filename.fits#1#   -- Specifies extension 1 for a multi-extension FITS file of binary tables

You can explore the data, perform simple analysis (e.g. EW measurements) and more. See the Notebook for more.

Warning

The xspec gui is still in an experimental stage and some aspects might not work as expected.

lt_plot

Plot one or more spectra. This has fewer features than lt_xspec above, but is faster.

For example:

lt_plot filename1 filename2

A plot of the spectrum in filename1 appears, and you can navigate around it using the same commands as in lt_xspec. Move to the next or previous spectrum with the right or left arrow key.

To list all the command line options available, use:

lt_plot -h

lt_absline

This plots a single absorption line, given a transition rest wavelength (Angstroms) or name (e.g. CIV1548), log10 column density, and Doppler parameter (km/s).

For example:

lt_absline 1215.6701 14.0 30
lt_absline HI1215 14.0 30

plots a Hydrogen Ly-a line with column density of 1014 cm-2 and b=30 km/s. A plot will appear and the line info and EW as well, i.e.

[AbsLine: HI 1215, wrest=1215.6700 Angstrom]
EW = 0.268851 Angstrom

Try:

lt_absline -h

for the full set of options.

lt_radec

Input coordinates in one format and print them out in several formats:

lt_radec 152.25900,7.22885
lt_radec J100902.16+071343.86
lt_radec 10:09:02.16,+07:13:43.8

lt_line

Print the atomic data for an input ion, transition or for an entire linelist.:

lt_line -h
lt_line HI
lt_line HI -z 3.0
lt_line HI1215
lt_line 1215
lt_line --all

lt_continuumfit

Launch the GUI to continuum fit a spectrum. If a redshift is supplied by zsys, then the script assumes this is a QSO.:

lt_continuumfit input_file output_filename --redshift 0.867

Here is the current usage message:

usage: lt_continuumfit [-h] [--redshift REDSHIFT] [--wchunk WCHUNK]
                       file outfil

GUI to fit a continuum to a spectrum

positional arguments:
  file                 Input spectral file (FITS, ASCII, etc.)
  outfil               Output, normalized spectrum filename; FITS [can be the
                       same]

optional arguments:
  -h, --help           show this help message and exit
  --redshift REDSHIFT  Redshift of the Source
  --wchunk WCHUNK      Width of a 'chunk' (Ang)

Reference & API

The linetools API

linetools.abund.ions Module

Utilities for working with ionized atoms.

Functions
ion_name(ion[, flg, nspace]) Convert ion tuple into a string
name_ion(ion) Convert string into ion tuple

linetools.abund.roman Module

Convert to and from Roman numerals.

This code was originally written by Mark Pilgrim (f8dy@diveintopython.org), and released under a Python 2.1.1 license, available at https://www.python.org/download/releases/2.1.1/license

Functions
fromRoman(s) Convert a Roman numeral to an integer.
toRoman(n) Convert an integer to Roman numeral.

linetools.abund.solar Module

Simple Solar abundance calculations.

Classes
SolarAbund([ref, verbose]) Class to handle simple Solar Abundance calculations

linetools.analysis.absline Module

Utlities for the analysis of absorption lines

Functions
aodm(spec, idata) AODM calculation on an absorption line
linear_clm(obj) Return N and sig_N given logN, sig_logN in an object
log_clm(obj) Return logN and sig_logN given linear N, sig_N
photo_cross(Z, ion, E[, datfil, silent]) Estimate photo-ionization cross-section using Fit parameters
sum_logN(obj1, obj2) Add log columns and return value and errors, with logic

linetools.analysis.cog Module

Curve of Growth calculations

Functions
cog_plot(COG_dict) Generate a plot for COG solution
single_cog_analysis(wrest, f, EW[, sig_EW, ...]) Perform COG analysis on a single component
Classes
single_cog_model Generate a single COG model
Class Inheritance Diagram

Inheritance diagram of linetools.analysis.cog.single_cog_model

linetools.analysis.continuum Module

Module for fitting a QSO continuum

Functions
Akima_co(wa, knots) Akima interpolation through the spline knots.
chisq_chunk(model, fl, er, masked, indices, ...) Calc chisq per chunk, update knots flags inplace if chisq is acceptable.
estimate_continuum(s, knots, indices, masked) Iterate to estimate the continuum.
find_continuum(spec[, edges, ax, debug, kind]) Estimate a continuum for a spectrum.
linear_co(wa, knots) linear interpolation through the spline knots.
make_chunks_qso(wa, redshift[, divmult, ...]) Generate a series of wavelength chunks for use by prepare_knots, assuming a QSO spectrum.
prepare_knots(wa, fl, er, edges[, ax, debug]) Make initial knots for the continuum estimation.
remove_bad_knots(knots, indices, masked, fl, er) Remove knots in chunks without any good pixels.
unmask(masked, indices, wa, fl, er[, minpix]) Forces each chunk to use at least minpix pixels.
update_knots(knots, indices, fl, masked) Calculate the y position of each knot.

linetools.analysis.interactive_plot Module

Classes for making interactive plots.

Functions
local_median(wa, fl, er, x[, npix, default]) find the median flux value at x using +/- npix pixels.
Classes
InteractiveCoFit(wa, fl, er, contpoints[, ...]) Class for interactively fitting a continuum
PlotWrapBase() A base class that has all the navigation and smoothing keypress events.
PlotWrapNav(fig, ax, wa, fl, artists[, ...]) Enable simple XIDL-style navigation for plotting a spectrum.
Class Inheritance Diagram

Inheritance diagram of linetools.analysis.interactive_plot.InteractiveCoFit, linetools.analysis.interactive_plot.PlotWrapBase, linetools.analysis.interactive_plot.PlotWrapNav

linetools.analysis.interp Module

Interpolation-related tools.

Functions
interp_Akima(x_new, x, y) Return interpolated data using Akima’s method.
Classes
AkimaSpline(xvals, yvals) Describes an Akima Spline through a set of points.

linetools.analysis.utils Module

Line analysis tools

These are intended to be methods generic to emission and absorption (e.g. Equivalent width)

Functions
box_ew(spec) Boxcar EW calculation
gaussian_ew(spec, ltype[, initial_guesses]) EW calculation using Gaussian fit

linetools.analysis.voigt Module

Voigt profiles

Heavily adapted from code by Ryan Cooke (e.g. alis)

Functions
voigt_from_abslines(iwave, line[, fwhm, ...]) Generates a Voigt model from a line or list of AbsLines
voigt_tau(wave, par) Find the optical depth at input wavelengths
voigt_wofz(vin, a) Uses scipy function for calculation.
voigtking(vin, a) Take from Alis by Ryan Cooke
Classes
single_voigt_model Generate a single Voigt model in astropy framework for fitting
Class Inheritance Diagram

Inheritance diagram of linetools.analysis.voigt.single_voigt_model

linetools.isgm.abscomponent Module

Classes for absorption line component

Classes
AbsComponent(radec, Zion, zcomp, vlim[, Ej, ...]) Class for an absorption component
Quantity A Quantity represents a number with some associated unit.
XSpectrum1D(wave, flux[, sig, co, units, ...]) A Class containing 1D spectra and associated methods

linetools.lists.linelist Module

Contains the LineList class

Classes
LineList(llst_key[, verbose, closest, ...]) This Class is designed to organize and handle information about atomic and/or molecular transition lines (e.g.
MaskedColumn Define a masked data column for use in a Table object.
Class Inheritance Diagram

Inheritance diagram of linetools.lists.linelist.LineList

linetools.lists.mk_sets Module

Used to make set lists. Only intended for developer use.

Functions
add_galaxy_lines(outfil[, infil, stop]) Append galaxy lines (as necessary)
mk_hi([infil, outfil, stop]) Make the HI list ISM + HI

linetools.lists.parse Module

Tools for parsing Line List data

Functions
grab_galaxy_linelists([do_this]) Pulls galaxy emission line lists from DESI project
line_data([nrows]) Defines the dict (and/or Table) for spectral line Data
mask_gal(data[, mask_keys]) Masks linelist attributes for all galaxy lines
mktab_morton00([do_this, outfil]) Used to generate a FITS Table for the Morton2000 paper
mktab_morton03([do_this, outfil, fits]) Used to generate a VO or FITS Table for the Morton2003 paper
open(filename[, mode, encoding, errors, ...]) Open an encoded file using the given mode and return a wrapped version providing transparent encoding/decoding.
parse_morton00([orig]) Parse tables from Morton 2000, ApJS, 130, 403
parse_morton03([orig, tab_fil, HIcombine]) Parse tables from Morton 2003, ApJS, 149, 205
parse_verner96([orig, write]) Parse tables from Verner, Verner, & Ferland (1996, Atomic Data and Nuclear Data Tables, Vol.
read_CO() Simple def to read CO UV data
read_H2() Simple def to read H2 data
read_euv() read additional EUV lines
read_forbidden() read galaxy emission lines (forbidden)
read_galabs() read galaxy absorption lines
read_recomb() read galaxy emission lines (recombination)
read_sets([infil]) Read sets file
read_verner94() Read Verner1994 Table
update_fval(table[, verbose]) Update f-values from the literature
update_gamma(table) Update/add-in gamma values
update_wrest(table[, verbose]) Update wrest values (and Ej,Ek)

linetools.spectra.convolve Module

Functions related to convolution.

Functions
convolve_psf(array, fwhm[, boundary, ...]) Convolve an array with a gaussian kernel.

linetools.spectra.io Module

Reading and writing of spectra

Functions
chk_for_gz(filenm) Checks for .gz extension to an input filename and returns file
get_cdelt_dcflag(hd) Find the wavelength stepsize and dcflag from a fits header.
get_table_column(tags, hdulist[, idx]) Find a column in a FITS binary table
get_wave_unit(tag, hdulist[, idx]) Attempt to pull wavelength unit from the Table
give_wv_units(wave) Give a wavelength array units of Angstroms, if unitless.
is_UVES_popler(hd) Check if this header is UVES_popler output.
parse_DESI_brick(hdulist[, select]) Read a spectrum from a DESI brick format HDU list
parse_FITS_binary_table(hdulist[, exten, ...]) Read a spectrum from a FITS binary table
parse_UVES_popler(hdulist) Read a spectrum from a UVES_popler-style fits file.
parse_hdf5(inp, **kwargs) Read a spectrum from HDF5 written in XSpectrum1D format
parse_linetools_spectrum_format(hdulist) Parse an old linetools-format spectrum from an hdulist
parse_two_file_format(specfil, hdulist[, efil]) Parse old two file format (one for flux, another for error).
readspec(specfil[, inflg, efil, verbose, ...]) Read a FITS file (or astropy Table or ASCII file) into a
setwave(hdr) Generate wavelength array from a header

linetools.spectra.lsf Module

Module for dealing with LSFs of various astronomical instruments.

Classes
LSF(instr_config) Class to deal with line-spread-functions (LSFs) from various different astronomical spectrographs.

linetools.spectra.plotting Module

Plotting utilities.

Functions
get_flux_plotrange(fl[, perc, mult_pos, ...]) Find y limits for a spectrum.

linetools.spectra.xspectrum1d Module

Module containing the XSpectrum1D class

Classes
Column Define a data column for use in a Table object.
Table([data, masked, names, dtype, meta, ...]) A class to represent tables of heterogeneous data.
UnitBase Abstract base class for units.
XSpectrum1D(wave, flux[, sig, co, units, ...]) A Class containing 1D spectra and associated methods
Class Inheritance Diagram

Inheritance diagram of linetools.spectra.xspectrum1d.XSpectrum1D

linetools.spectralline Module

Classes for an emission or absorption spectral line

Functions
many_abslines(all_wrest, llist) Generate a list of AbsLine objects.
Classes
AbsLine(trans, **kwargs) Class representing a spectral absorption line.
EmLine(trans, **kwargs) Class representing a spectral emission line.
LineLimits(wrest, z, zlim, **kwargs) An object for handling the ‘limits’ of a line
SpectralLine(ltype, trans[, linelist, ...]) Class for a spectral line.
XSpectrum1D(wave, flux[, sig, co, units, ...]) A Class containing 1D spectra and associated methods
Class Inheritance Diagram

Inheritance diagram of linetools.spectralline.AbsLine, linetools.spectralline.EmLine, linetools.spectralline.SpectralLine

linetools.utils Module

Module for general utilities which don’t belong in another sub-package.

Functions
between(a, vmin, vmax) Return a boolean array True where vmin <= a < vmax.
convert_quantity_in_dict(idict) Return a dict where Quantities (usually from a JSON file)
give_dv(z, zmean[, rel]) Gives velocity difference between z and zmean, at zmean.
give_dz(dv, zmean[, rel]) Gives redshift difference for a given velocity difference(s) at zmean.
jsonify(obj[, debug]) Recursively process an object so it can be serialised in json format.
loadjson(filename) Load a python object saved with savejson.
overlapping_chunks(chunk1, chunk2) True if there is overlap between chunks chunk1 and chunk2.
radec_to_coord(radec) Converts one of many of Celestial Coordinates radec formats to an astropy SkyCoord object.
rel_vel(wavelength, wv_obs) Simple relative velocity method
savejson(filename, obj[, overwrite, indent, ...]) Save a python object to filename using the JSON encoder.
scipy_rebin(a, *args) Simple script to rebin an input array to a new shape.
v_from_z(z1, z2) Find the relativistic velocity between 2 redshifts.
z_from_v(z, v) Find the redshift given z and v

Indices and tables