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.11 or later
`astropy`_ version 1.3 or later
scipy version 0.16 or later
matplotlib version 1.4 or later
PyQT5 version 5 (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 h5py

For PyQT, the current conda version is PyQT5. If you are still using PyQT4, then you will need to update to use linetools GUIs:

conda install pyqt=5

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¶

There is currently a pip wheel on PyPi but it is woefully out of date. We will try to update before long. But for now follow the instructions in the section below, Installing Linetools from Source to install linetools.

Also note, 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.3 (Unreleased)¶

Updates¶

Added extra attributes to AbsComponents
Added script to get HST/COS life-time position from date
Added Cashman+2017 catalog in LineList
Significant refactor of AbsComponent
LineList “AGN” added

Bug fixes¶

0.2 (2017-10-24)¶

Updates¶

Extra features to main Objects like XSpectrum1D, AbsComponent, AbsLine, LineList
Added some extra emisison lines to Galaxy LineList
Refactor from pyQt4 -> pyQt5
Improvements to GUIs and scripts
Added EmLine and EmSystem classes
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 XSpectrum1D.get_local_s2n, a method to calculate average signal-to-noise at a given wavelength
Added xabssysgui GUI
Added new LineList(e.g. Galaxy)
Added EmLine child to SpectralLine
Added LineLimits class
Added SolarAbund class
Added lt_radec and lt_line scripts
Added LSF class to handle line-spread-functions. Currently implemented for HST/COS and most HST/STIS configurations.

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

wrest	z	flag_N	logN	sig_logN
Angstrom
float64	float64	int64	float64	float64
1215.67	2.92939	1	14.1	0.15
1025.7222	2.92939	1	14.15	0.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

Z	ion	A	Ej	z	vmin	vmax	flag_N	logN	sig_logN
					km / s	km / s
int64	int64	int64	float64	float64	float64	float64	int64	float64	float64
1	1	0	0.0	2.92939	-300.0	300.0	1	14.1172024817	0.117911610801
14	2	0	0.0	2.92939	-300.0	300.0	1	13.9006157733	0.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()

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'])

4221 Angstrom 1 13.5883729709 0.0150745701489
3702 Angstrom 1 13.7708705955 0.0862006463782
707 Angstrom 1 13.6707360009 0.0640855113383
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	zcomp	float	Absorption redshift of the component
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 ~linetools.spectralline.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.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 (e.g. a JSON file):

abscomp = AbsComponent.from_dict(idict)

One may also generate a set of components from a larger list of AbsLines:

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

Attributes¶

There is a set of default attributes that are initialized at Instantiation and kepy in an attrib dict. These are:

init_attrib = {'N': 0./u.cm**2, 'sig_N': [0.,0.]/u.cm**2, 'flag_N': 0, # Column
          'logN': 0., 'sig_logN': np.array([0.,0.]),
          'b': 0.*u.km/u.s, 'sig_b': 0.*u.km/u.s,  # Doppler
          'vel': 0*u.km/u.s, 'sig_vel': 0*u.km/u.s
          }

One can access these attributes with standard . syntax, e.g.:

logN = abscomp.logN   # This accesses the logN value from the internal attrib dict

Inspecting¶

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

Generate a Table¶

If the class contains one or more AbsLines, you may generate a ~astropy.table.Table 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 column densities are synthesized at instantiation unless one sets skip_synth=True.

There are two approaches to synthesis:

Standard¶

The default uses the synthesize_colm() method. Positive, unsaturated detections are combined in a weighted mean whereas limits are compared and the most sensitive one is adopted.:

Here is the set of rules:

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

A future implementation will introduce a flag for bracketing an upper and lower limit value.

Median¶

Another option for synthesizing the column densities and other attributes is to take the median values. This is performed when one passes adopt_median=True to the from_abslines() call, e.g.:

abscomp = AbsComponent.from_abslines([lya,lyb], adopt_median=True)

The current code computes the median on all input values and does not consider the flag_N values for the input AbsLine objects.

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 AbsComponent with the to_dict() method:

cdict = component.to_dict()

One may also wish to Voigt profile fit components with one of a number of software packages (e.g. ALIS, JoeBVP). To generate an input file for JoeBVP use:

from linetools.isgm.io import write_joebvp_from_components
write_joebvp_from_components(component_list, specfile, 'output_file.ascii')

Similarly, one can generate a list of components from an outputted JoeBVP file:

from linetools.isgm.io import read_joebvp_to_components
comp_list = read_joebvp_to_components('joebvp_file.out')

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.

Tables¶

You can also create a list of components using an input astropy.table.Table object (mandatory column names are [‘RA’, ‘DEC’, ‘ion_name’, ‘z_comp’, ‘vmin’, ‘vmax’, ‘Z’, ‘ion’, ‘Ej’]):

from astropy.table import Table
tab = Table()
tab['ion_name'] = ['HI', 'HI', 'CIV', 'SiII', 'OVI']
tab['Z'] = [1,1,4,14,8]
tab['ion'] = [1,1,4,2,6]
tab['Ej'] = [0.,0.,0.,0.,0.]/u.cm
tab['z_comp'] = [0.05, 0.0999, 0.1, 0.1001, 0.3]
tab['logN'] = [12.0, 13.0., 13.0, 14.0, 13.5]
tab['sig_logN'] = [0.1, 0.12., 0.15, 0.2, 0.1]
tab['flag_logN'] = [1, 1, 2, 0, 1]
tab['RA'] = 100.0 * u.deg * np.ones(len(tab))
tab['DEC'] = -0.8 * u.deg * np.ones(len(tab))
tab['vmin'] = -50 * u.km/u.s * np.ones(len(tab))
tab['vmax'] =  50 * u.km/u.s * np.ones(len(tab))

complist = ltiu.complist_from_table(tab)

These components will not have AbsLines defined by default. However, it is very easy to append AbsLines to these components for a given observed wavelength and minimum rest-frame equivalent width:

# add abslines
wvlim = [1100, 1500]*u.AA
for comp in complist:
   comp.add_abslines_from_linelist(llist='ISM', wvlim=wvlim, min_Wr=0.01*u.AA, chk_sep=False)

This will look for transitions in the LineList(‘ISM’) object and append those expected to be within wvlim to each corresponding AbsComponent in the complist. In this example, only those AbsLines expected to have rest-frame equivalent widths larger than 0.01*u.AA will be appended, but if you wish to include all of the available AbsLines you can set min_Wr=None. Here, we have set chk_sep=False to avoid checking for coordinates because by construction the individual AbsLines have the same coordinates as the corresponding AbsComponent.

One can, of course, go the other way and generate a Table from a list of components:

tab2 = ltiu.table_from_complist(complist)

If you wish to write to disk, we recommend that you use the astropy format: ascii.ecsv

Use Components to create a spectrum model¶

You can create a spectrum model using a list of AbsComponents, e.g.:

from linetools.analysis import voigt as lav
wv_array = np.array(1100,1500, 0.1)*u.AA
model = lav.voigt_from_components(wv_array, complist)

In this manner, model is a XSpectrum1D object with the the AbsComponents contained in the complist list. You may also add a convolution with a given kernel to compare with observations from different spectrographs, in which case you can use:

model = lav.voigt_from_components(wv_array, complist, fwhm=3)

This is same as before but the model is convolved with a Gaussian kernel of FWHM of 3 pixels.

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.

List¶

One can generate a list of all the AbsLine objects with:

lines = abssys.list_of_abslines()

get absline¶

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

measure_restew¶

Measure the rest-frame equivalent widths for all AbsLine objects in the system. In addition to the inputs described in measure_ew, each line to be measured must have line.analy[‘do_anlaysis’]=1. Here is an example:

abssys.measure_restew(spec=xspec_object)

fill transitions¶

Generate an astropy.Table of information on the absorption lines in the system. This is stored in the ._trans attribute:

abssys.fill_trans()
print(abssys._trans)

Components¶

~linetools.igsm.abssystem.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))

~linetools.igsm.abssystem.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.

AbsSightline Class¶

Notebooks¶

Overview¶

This Class is designed to organize the absorption systems along a single sightline. It may be most commonly used for extragalactic sightlines, but it can be applied to the ISM as well.

By definition, an AbsSightline is a unique collection of absorption components. The only quantities required

to define the AbsSightline are its coordinates on the sky.

Instantiation¶

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

abssl = GenericAbsSightline((10.0*u.deg, 45*u.deg))

More commonly, one will instantiate with one a set of components:

lya = AbsLine('HI 1215', z=2.3)
 lya.limits.set([-300.,300.]*u.km/u.s)  # vlim
 lyb = AbsLine(1025.7222*u.AA, z=2.3)
 lyb.limits.set([-300.,300.]*u.km/u.s)  # vlim
 abscomp = AbsComponent.from_abslines([lya,lyb])
 abscomp.coord = ltu.radec_to_coord((10.*u.deg, 45*u.deg))
 abssl = GenericAbsSightline.from_components([abscomp])

Inspecting¶

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

Generate a Table¶

If the class contains one or more AbsComponent obejcts, you may generate a ~astropy.table.Table from their attributes and data:

comp_tbl = abssl.build_table()

I/O¶

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

asldict = abssl.to_dict()

This can then be written to disk with a JSON or yaml dump.

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:

>>> print(sol['C'])
8.43
>>> print(sol[6])
8.43

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

Multiple elements can also be selected:

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

Element ratios can be accessed using the get_ratio method:

>>> print(sol.get_ratio('C/Fe')) 
0.98

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
‘AGN’: Lines tipically identified in AGN 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 gives 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 Table 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>}

Alternatively, we could call this method via tuple, e.g., (8,6), where the first entry is the atomic number (of O) and the second is the ionization state (VI)

ovi = ism.all_transitions((8,6))
print(ovi['name', 'wrest', 'f'])
  name     wrest     f
          Angstrom
-------- --------- ------
OVI 1031 1031.9261 0.1325
OVI 1037 1037.6167 0.0658

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', line=889.7228*u.AA, or line=(8,6). 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
Limits	limits	LineLimits	The redshift and limits of the line in redshift, velocity (w/r to its redshift) and observed wavelength.
RA, Dec	attrib[‘coord’]	Coord	astropy.coordinate
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

Note that redshift is sufficiently important that it is contained within its own object. It is also accessible as a property, e.g.:

z = sline.z

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 above) and the LineLimits of the line must have previously been specified.

When executed, the EW and sig_EW attibutes are filled:

specline.measure_ew()
print(specline.attrib['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 other 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¶

~linetools.spectra.xspectrum1d.XSpectrum1D contains and manipulates 1-d spectra, each of 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 at instantiation. 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 ~linetools.spectra.xspectrum1d.XSpectrum1D.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 # doctest: +SKIP
<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 spectra from specdb), 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 ~linetools.spectra.io.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 (masking=’None’) 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¶

~linetools.spectra.xspectrum1d.XSpectrum1D.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¶

~linetools.spectra.xspectrum1d.XSpectrum1D.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 ~linetools.spectra.xspectrum1d.XSpectrum1D.normalize method. This also sets spec.normed to True.

Finally, you can apply small variations to the continuum anchor points with ~linetools.spectra.xspectrum1d.XSpectrum1D.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 ~linetools.spectra.xspectrum1d.XSpectrum1D.box_smooth, ~linetools.spectra.xspectrum1d.XSpectrum1D.gauss_smooth, and ~linetools.spectra.xspectrum1d.XSpectrum1D.ivar_smooth.

Other methods¶

You can join one XSpectrum1D instance with another overlapping spectrum using ~linetools.spectra.xspectrum1d.XSpectrum1D.splice. ~linetools.spectra.xspectrum1d.XSpectrum1D.pix_minmax finds the pixel indices corresponding to a wavelength or velocity range, and ~linetools.spectra.xspectrum1d.XSpectrum1D.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 (~linetools.spectra.xspectrum1d.XSpectrum1D.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: ~linetools.spectra.xspectrum1d.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 ~linetools.spectra.io.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¶

xspec Documentation¶

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")

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 the command line (recommended)¶

We recommend you use the script provided with linetools.

Then it is as simple as:

> lt_xspec filename

Here are the current command-line options:

> lt_xspec -h
usage: lt_xspec [-h] [--zsys ZSYS] [--norm] [--exten EXTEN]
                [--wave_tag WAVE_TAG] [--flux_tag FLUX_TAG]
                [--sig_tag SIG_TAG] [--var_tag VAR_TAG] [--ivar_tag IVAR_TAG]
                file

Parse for XSpec

positional arguments:
  file                 Spectral file

optional arguments:
  -h, --help           show this help message and exit
  --zsys ZSYS          System Redshift
  --norm               Show spectrum continuum normalized (if continuum is
                       provided)
  --exten EXTEN        FITS extension
  --wave_tag WAVE_TAG  Tag for wave in Table
  --flux_tag FLUX_TAG  Tag for flux in Table
  --sig_tag SIG_TAG    Tag for sig in Table
  --var_tag VAR_TAG    Tag for var in Table
  --ivar_tag IVAR_TAG  Tag for ivar in Table

From within ipython or equivalent¶

from linetools.guis import xspecgui as ltxsg

import imp; imp.reload(ltxsg)
ltxsg.main(spec_fil)

Navigating - These key strokes help you explore the spectrum (be sure to click in the spectrum panel first!)¶

Setting the window edges (mouse+keystroke)¶

l – Set left edge of window
r – Set right edge of window
t – Set top edge of window
b – Set bottom edge of window
Z – Set bottom edge to 0.
W – View full spectrum
s,s – Set a zoom-in window at 2 mouse positions

Zoom in/out Wavelength¶

i – Zoom in on cursor
I – Zoom in extra fast
o – Zoom out
O – Zoom out extra fast

Best estimate Flux limits¶

y – Make a guess for the flux plotting limits

Zoom out Flux¶

Y – Zoom out

Pan¶

[ – Pan left
{ – Pan left extra
] – Pan right
} – Pan right extra

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 the command line (recommended)¶

We recommend you use the script provided with linetools.

Then it is as simple as:

> lt_xabssys spec_file abssys_file

The abssys_file is expected to be a JSON file that contains an AbsSystem (likely written with the write_json method).

Here are the current command-line options:

> lt_xspec -h
usage: lt_xabssys [-h] [-outfile OUTFILE] [-llist LLIST] [--un_norm]
              spec_file abssys_file

Parse for XAbsSys

positional arguments:
  spec_file         Spectral file
  abssys_file       AbsSys file (JSON)

optional arguments:
  -h, --help        show this help message and exit
  -outfile OUTFILE  Output filename
  -llist LLIST      Name of LineList
  --un_norm         Spectrum is NOT normalized

From within ipython or equivalent¶

Not yet implemented

Navigating in the Main Window- These key strokes help you explore the spectrum (be sure to click in the spectrum panel first!)¶

Setting the window edges (mouse+keystroke)¶

l – Set left edge of window
r – Set right edge of window
t – Set top edge of window
b – Set bottom edge of window

Zoom in/out Wavelength¶

i – Zoom in on cursor
I – Zoom in extra fast
o – Zoom out
O – Zoom out extra fast

Zoom out Flux¶

Y – Zoom out

Pan¶

[ – Pan left
{ – Pan left extra
] – Pan right
} – Pan right extra

Adjusting Rows/Columns¶

= – Move to next page of lines
- – Move to previous page of lines
C – Add a column
c – Remove a column
K – Add a row
k – Remove a row

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.

Here is the full usage:

usage: lt_xspec [-h] [-z ZSYS] [--norm] [--air] [--exten EXTEN]
                [--splice SPLICE] [--scale SCALE]
                file

Parser for lt_xspec v1.2; Note: Extra arguments are passed to read_spec (e.g.
--flux_tag=FX)

positional arguments:
  file                  Spectral file; specify extension by appending #exten#

optional arguments:
  -h, --help            show this help message and exit
  -z ZSYS, --zsys ZSYS  System Redshift
  --norm                Show spectrum continuum normalized (if continuum is
                        provided)
  --air                 Convert input spectrum wavelengths from air to vacuum
  --exten EXTEN         FITS extension
  --splice SPLICE       Splice with the input file; extension convention
                        applies
  --scale SCALE         Scale factor for GUI size [1. is default]

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 10¹⁴ 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

Here is the usage:

beast> lt_line -h
usage: lt_line [-h] [-a] [--llist LLIST] [-t TOLER] [-z REDSHIFT] [inp]

Print spectral line data of a line or lines.

positional arguments:
  inp                   Ion, transition name, or Rest wavelength (e.g. HI,
                        HI1215 or 1215)

optional arguments:
  -h, --help            show this help message and exit
  -a, --all             Print all lines
  --llist LLIST         Name of LineList: ISM, HI, H2, CO, etc.
  -t TOLER, --toler TOLER
                        Matching tolerance (in Ang) on input wavelength
  -z REDSHIFT, --redshift REDSHIFT
                        Matching tolerance (in Ang) on input wavelength

lt_solabnd¶

Print the solar abundances by number relative to Hydrogen on the traditional log scale of 12, e.g. e(Fe) = 7.55:

lt_solabnd Fe
lt_solabnd -a
lt_solabnd -a --sortZ

Here is the usage:

beast> lt_solabnd -h
usage: lt_solabnd [-h] [-a] [--sortZ] [inp]

Print Solar abundance data for an element or all elements.

positional arguments:
  inp         Elm (e.g. H, Fe)

optional arguments:
  -h, --help  show this help message and exit
  -a, --all   Print all values
  --sortZ     Sort on Atomic Number

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)