I am preparing for some spectral analysis which will involve combining three overlapping spectra. Two of the spectra have "dimensionless" units (one, and albedo, has units of flux/flux; the other, a filter response, is in photons/photons). I would like to use specutils and astropy to make this easier and generally avoid reinventing any wheels.
However, I am having difficulty converting my data into a Spectrum1D-type object. Here is the relevant code snippet:
import numpy as np
import astropy
from astropy import units as u
from specutils import Spectrum1D
dataSpectrumFileName = 'filename.dat' # Declare filename
dataSpectrumRaw = list(np.loadtxt(dataSpectrumFileName).T) # Load file.
#Now dataSpectrumRaw[0] is the wavelength grid and dataSpectrumRaw[1] is the dimensionless spectrum.
# Declare units
dataSpectrumRaw[0] = [ u.Quantity(dataSpectrumRaw[0][i],u.micrometer) for i in
range(len(dataSpectrumRaw[0])) ]
dataSpectrumRaw[1] = [ u.Quantity(dataSpectrumRaw[1][i],u.dimensionless_unscaled) for i in
range(len(dataSpectrumRaw[1])) ]
# Convert to Spectrum1D item
dataSpectrum = Spectrum1D(spectral_axis = dataSpectrumRaw[0], flux = dataSpectrumRaw[1])
My expectation is that this will load in my data, convert the raw numerical information into two lists of Quantity-type objects (the former with units of u.micrometer and the latter with units of u.dimensionless_unscaled), then produce a Spectrum1D-type object named dataSpectrum with all the relevant information and associated functionality. Indeed, a check prior to running the final line of my code snippet reveals:
>>>dataSpectrumRaw[0][0]
0.2 um
>>>type(dataSpectrumRaw[0][0])
<class 'astropy.units.quantity.Quantity'>
>>>dataSpectrumRaw[1][0]
0.269681
>>>type(dataSpectrumRaw[1][0])
<class 'astropy.units.quantity.Quantity'>
However, on compiling the above code entirely I am given an error:
Traceback (most recent call last):
File "/Users/[REDACTED]/specutilsExperiments.py", line 40, in <module>
dataSpectrum = Spectrum1D(spectral_axis = dataSpectrumRaw[0], flux = dataSpectrumRaw[1])
File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/specutils/spectra/spectrum1d.py", line 71, in __init__
raise ValueError("Flux must be a `Quantity` object.")
ValueError: Flux must be a `Quantity` object.
The documentation and error text both seem to indicate that what I am trying to do SHOULD work. I have also tried converting whole arrays into u.Quantity objects, which for some reason does not stick (the elements convert back to np.float64 before I can operate on them further).
I have then two questions: 1) Why does this not work as expected? 2) How can I make it work as expected?
Please answer both questions if possible; without the answer to 1) I will likely run into this again.
I think you're making things more complicated for yourself than necessary. Try it like this:
>>> spectrum_raw = np.loadtxt(filename).T # no reason to convert to a list
>>> spectral_axis = spectrum_raw[0] * u.micrometer
>>> flux = spectrum_raw[1] * u.dimensionless_unscaled
>>> spectrum = Spectrum1D(spectral_axis=spectral_axis, flux=flux)
>>> spectrum
<Spectrum1D(flux=<Quantity [0.79502858, 0.68160812, 0.04951287, 0.86626975, 0.35146137,
0.13505172, 0.24077603, 0.65288882, 0.24230813, 0.51846452]>, spectral_axis=<Quantity [0.41284773, 0.23807243, 0.14288542, 0.99212095, 0.75582049,
0.1175328 , 0.66007513, 0.48731926, 0.07727119, 0.02389908] um>)>
In brief, you are creating Python lists of Quantity scalars, where list is the standard heterogeneous collection in Python with no specific meaning. But an Astropy Quantity can store either a single scalar value, or an arbitrarily-sized array of values of the same unit. So that's what you want when you already have an array of values in the same units. See Creating Quantity instances.
It seems maybe you tried this at some point:
I have also tried converting whole arrays into
u.Quantityobjects, which for some reason does not stick.
Without knowing what you did I can't say why it did not "stick", but the above example should work. Not also that you can convert an ndarray to a Quantity simply my multiplying the array by the desired unit. It is basically syntactic sugar equivalent to Quantity(array, unit).