1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
|
###########################
Working with Spectral Cubes
###########################
Spectral cubes can be read directly with :class:`~specutils.Spectrum`.
A specific example of this is demonstrated here. In addition to the functions
demonstrated below, as a subclass of `NDCube <https://github.com/sunpy/ndcube>`_,
:class:`~specutils.Spectrum` also inherits useful methods for e.g. cropping
based on combinations of world and spectral coordinates. Most of the
functionality inherited from `~ndcube.NDCube` requires initializing the
``Spectrum`` object with a WCS describing the coordinates for all axes of
the data.
Note that the workflow described here is for spectral cubes that are rectified
such that one of the axes is entirely spectral and all the spaxels have the same
``spectral_axis`` values (i.e., case 2 in :ref:`specutils-representation-overview`).
For less-rectified cubes, pre-processing steps (not addressed by specutils at the
time of this writing) will be necessary to rectify the cubes into that form.
Note, though, that such cubes can be stored in specutils data structures (cases
3 and 4 in :ref:`specutils-representation-overview`), which support *some* of
these behaviors, even though the fulkl set of tools do not yet apply.
Loading a cube
==============
We'll use a ``MaNGA cube`` for our example, and load the data from the
repository directly into a new ``Spectrum`` object:
.. code-block:: python
>>> from astropy.utils.data import download_file
>>> from specutils.spectra import Spectrum
>>> filename = "https://stsci.box.com/shared/static/28a88k1qfipo4yxc4p4d40v4axtlal8y.fits"
>>> file = download_file(filename, cache=True) # doctest: +REMOTE_DATA
>>> sc = Spectrum.read(file, format='MaNGA cube') # doctest: +REMOTE_DATA
The cube has 74x74 spaxels with 4563 spectral axis points in each one:
.. code-block:: python
>>> sc.shape # doctest: +REMOTE_DATA
(4563, 74, 74)
Print the contents of 3 spectral axis points in a 3x3 spaxel array:
.. code-block:: python
>>> sc[2000:2003, 30:33, 30:33] # doctest: +REMOTE_DATA
<Spectrum(flux=[[[0.4892023205757141 ... 0.5994223356246948]]] 1e-17 erg / (Angstrom s spaxel cm2) (shape=(3, 3, 3), mean=0.54165 1e-17 erg / (Angstrom s spaxel cm2)); spectral_axis=<SpectralAxis
(observer to target:
radial_velocity=0.0 km / s
redshift=0.0)
[5.73984286e-07 5.74116466e-07 5.74248676e-07] m> (length=3); uncertainty=InverseVariance)>
Spectral slab extraction
========================
The `~specutils.manipulation.spectral_slab` function can be used to extract
spectral regions from the cube.
.. code-block:: python
>>> import astropy.units as u
>>> from specutils.manipulation import spectral_slab
>>> ss = spectral_slab(sc, 5000.*u.AA, 5003.*u.AA) # doctest: +REMOTE_DATA
>>> ss.shape # doctest: +REMOTE_DATA
(3, 74, 74)
>>> ss[::, 30:33,30:33] # doctest: +REMOTE_DATA
<Spectrum(flux=[[[0.6103081107139587 ... 0.936118483543396]]] 1e-17 erg / (Angstrom s spaxel cm2) (shape=(3, 3, 3), mean=0.83004 1e-17 erg / (Angstrom s spaxel cm2)); spectral_axis=<SpectralAxis
(observer to target:
radial_velocity=0.0 km / s
redshift=0.0)
[5.00034537e-07 5.00149688e-07 5.00264865e-07] m> (length=3); uncertainty=InverseVariance)>
Spectral Bounding Region
========================
The `~specutils.manipulation.extract_bounding_spectral_region` function can be used to
extract the bounding region that encompases a set of disjoint `~specutils.SpectralRegion`
instances, or a composite instance of `~specutils.SpectralRegion` that contains
disjoint sub-regions.
.. code-block:: python
>>> from specutils import SpectralRegion
>>> from specutils.manipulation import extract_bounding_spectral_region
>>> composite_region = SpectralRegion([(5000*u.AA, 5002*u.AA), (5006*u.AA, 5008.*u.AA)])
>>> sub_spectrum = extract_bounding_spectral_region(sc, composite_region) # doctest: +REMOTE_DATA
>>> sub_spectrum.spectral_axis # doctest: +REMOTE_DATA +FLOAT_CMP
<SpectralAxis
(observer to target:
radial_velocity=0.0 km / s
redshift=0.0)
[5.00034537e-07, 5.00149688e-07, 5.00264865e-07, 5.00380068e-07,
5.00495298e-07, 5.00610555e-07, 5.00725838e-07] m>
Moments
=======
The `~specutils.analysis.moment` function can be used to compute moments of any order
along one of the cube's axes. By default, ``axis='spectral'``, in which case the moment
is computed along the spectral axis.
.. code-block:: python
>>> from specutils.analysis import moment
>>> m = moment(sc, order=1) # doctest: +REMOTE_DATA
>>> m.shape # doctest: +REMOTE_DATA
(74, 74)
>>> m[30:33,30:33] # doctest: +REMOTE_DATA +FLOAT_CMP
<Quantity [[6.97933331e-07, 6.98926463e-07, 7.00540974e-07],
[6.98959625e-07, 7.00280655e-07, 7.03511823e-07],
[7.00740294e-07, 7.04527986e-07, 7.08245958e-07]] m>
Use Case
========
Example of computing moment maps for specific wavelength ranges in a
cube, using `~specutils.manipulation.spectral_slab` and
`~specutils.analysis.moment`.
.. plot::
:include-source:
:align: center
:context: close-figs
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.utils.data import download_file
from specutils import Spectrum, SpectralRegion
from specutils.analysis import moment
from specutils.manipulation import spectral_slab
filename = "https://stsci.box.com/shared/static/28a88k1qfipo4yxc4p4d40v4axtlal8y.fits"
fn = download_file(filename, cache=True)
spec1d = Spectrum.read(fn)
# Extract H-alpha sub-cube for moment maps using spectral_slab
subspec = spectral_slab(spec1d, 6745.*u.AA, 6765*u.AA)
ha_wave = subspec.spectral_axis
# Extract wider sub-cube covering H-alpha and [N II] using spectral_slab
subspec_wide = spectral_slab(spec1d, 6705.*u.AA, 6805*u.AA)
ha_wave_wide= subspec_wide.spectral_axis
# Convert flux density to microJy and correct negative flux offset for
# this particular dataset
ha_flux = (np.sum(subspec.flux.value, axis=(1,2)) + 0.0093) * 1.0E-6*u.Jy
ha_flux_wide = (np.sum(subspec_wide.flux.value, axis=(1,2)) + 0.0093) * 1.0E-6*u.Jy
# Compute moment maps for H-alpha line
moment0_halpha = moment(subspec, order=0)
moment1_halpha = moment(subspec, order=1)
# Convert moment1 from AA to velocity
# H-alpha is redshifted to 6755 AA for this galaxy
print(moment1_halpha[40,40])
vel_map = 3.0E5 * (moment1_halpha.value - 6.755E-7) / 6.755E-7
# Plot results in 3 panels (subspec_wide, H-alpha line flux, H-alpha velocity map)
f,(ax1,ax2,ax3) = plt.subplots(1, 3, figsize=(15, 5))
ax1.plot(ha_wave_wide, (ha_flux_wide)*1000.)
ax1.set_xlabel('Angstrom', fontsize=14)
ax1.set_ylabel('uJy', fontsize=14)
ax1.tick_params(axis="both", which='major', labelsize=14, length=8, width=2, direction='in', top=True, right=True)
ax2.imshow(moment0_halpha.value, origin='lower')
ax2.set_title('moment = 0')
ax2.set_xlabel('x pixels', fontsize=14)
ax3.imshow(vel_map, vmin=-100., vmax=100., cmap='rainbow', origin='lower')
ax3.set_title('moment = 1')
ax3.set_xlabel('x pixels', fontsize=14)
|
