The goal of this notebook is to demonstrate how to read in and manipulate DESI spectra using simulated spectra created as part of a DESI Data Challenge.
If you identify any errors or have requests for additional functionality please create a new issue on https://github.com/desihub/desispec/issues or send a note to desi-data@desi.lbl.gov.
We'll assume that you've gone through the installation instructions at:
https://desi.lbl.gov/trac/wiki/Pipeline/GettingStarted/Laptop/JuneMeeting
as far as Try one of our tutorials.
We'll also assume that you've grabbed the DESI data team's example spectral files and put them in the directory
$DESI_ROOT/spectro/redux/dc17a2
If you don't have these yet, a thumb drive is being passed around, but they're also at NERSC in the directory
/scratch2/scratchdirs/sjbailey/desi/dc17a/spectro/redux/dc17a2
If you've checked out desispec
from git, this jupyter notebook lives in desispec/doc/nb
and you can execute it as jupyter notebook Intro_to_DESI_spectra.ipynb
Now, let's import all the modules we'll need:
In [1]:
import os
import numpy as np
import healpy as hp
from glob import glob
import fitsio
from collections import defaultdict
from desitarget import desi_mask
import matplotlib.pyplot as plt
%pylab inline
If any of these imports fail, you should go back through the installation instructions as far as Try one of our tutorials.
In [2]:
%set_env DESI_SPECTRO_REDUX=/Users/dummy/desi/spectro/redux/
%set_env SPECPROD=dc17a2
In [3]:
import os
def check_env():
for env in ('DESI_SPECTRO_REDUX', 'SPECPROD'):
if env in os.environ:
print('{} environment set to {}'.format(env, os.getenv(env)))
else:
print('Required environment variable {} not set!'.format(env))
In [4]:
check_env()
In [5]:
basedir = os.path.join(os.getenv("DESI_SPECTRO_REDUX"),os.getenv("SPECPROD"),"spectra-64")
subdir = os.listdir(basedir)
print(basedir)
print(subdir)
In [6]:
basedir = os.path.join(basedir,subdir[0])
subdir = os.listdir(basedir)
pixnums = np.array([int(pixnum) for pixnum in subdir])
print(basedir)
print(subdir)
In [7]:
basedir = os.path.join(basedir,subdir[0])
subdir = os.listdir(basedir)
print(basedir)
print(subdir)
The directory structure is:
$DESI_SPECTRO_REDUX/$SPECTRO/spectra-{nside}/{group}/{pix}/*-{nside}-{pix}.fits`
where group = nside//100
. For example for nside=64
and pixel=16879
:
$DESI_SPECTRO_REDUX/$SPECTRO/spectra-64/168/16879/spectra-64-16879.fits
$DESI_SPECTRO_REDUX/$SPECTRO/spectra-64/168/16879/zbest-64-16879.fits
where the first file contains the spectra and the second file contains information on the best-fit redshifts from the redrock code.
In this directory-and-file structure, nside
is the HEALPix resolution and pixel
is the pixel number at that resolution. If you aren't familiar with HEALPix, it is an equal-area splitting of the sphere, where the sphere is initially divided into 12 equal-area pixels, and then each of those pixels is divided into 4 new equal-area pixels as nside
increases (a quad tree). Schematically, here's how nside
corresponds to pixel area in degrees:
In [8]:
sphere_area = 4*180.*180./np.pi
hpx_area = sphere_area
for i in range(10):
nside = 2**i
if i == 0:
hpx_area/=12.
else:
hpx_area/=4.
print(nside,hpx_area)
The nside
at which the example spectra are grouped therefore corresponds to ~0.84 sq. deg. Note that I could have checked this more easily (but less pedagogically) using the useful python HEALPix library:
In [9]:
hp.nside2pixarea(64,degrees=True)
Out[9]:
The spectra are stored in this fashion so that they are grouped (roughly) contiguously on the sky, with a reasonable number of spectra in each directory. It's easy to derive the approximate RA/Dec near each pixel number (note that we sneakily stored the pixel numbers as pixnums
when we were examining the directory structure):
In [10]:
ras, decs = hp.pix2ang(64, pixnums, nest=True, lonlat=True)
Note that the DESI Data Model will always use the NESTED scheme for HEALPix.
In [11]:
zipper = zip(pixnums,ras,decs)
[print("Pixel(nside=64): {} RA: {} DEC: {}".format(pix,ra,dec)) for pix,ra,dec in zipper]
Out[11]:
What about the Data Model for the spectra themselves? Let's take a look (remember that we built a directory containing a spectrum as basedir
above):
In [12]:
specfilename = glob(os.path.join(basedir,'spectra*fits'))[0]
DM = fitsio.FITS(specfilename)
DM
Out[12]:
The first extension FIBERMAP
stores the mapping of the imaging information used to target and place a fiber on the source:
In [13]:
fm = fitsio.read(specfilename,1)
fm.dtype.descr
Out[13]:
TARGETID
is the unique mapping from target information to a fiber. So, if you wanted to look up full imaging information for a spectrum, you can map back to target files using TARGETID
.
Just out of interest, are the RAs and Decs of these objects in the expected HEALPix pixel?
In [14]:
pixnums = hp.ang2pix(64, fm["RA_TARGET"], fm["DEC_TARGET"], nest=True, lonlat=True)
print(np.min(pixnums),np.max(pixnums))
print(specfilename)
In [15]:
plt.plot(fm["RA_TARGET"],fm["DEC_TARGET"],'b.')
Out[15]:
In [16]:
DM
Out[16]:
The remaining extensions store the wavelength, flux, inverse variance on the flux, mask and resolution matrix for the B, R and Z arms of the spectrograph. Let's determine the wavelength coverage of each spectrograph:
In [17]:
bwave = fitsio.read(specfilename,"B_WAVELENGTH")
rwave = fitsio.read(specfilename,"R_WAVELENGTH")
zwave = fitsio.read(specfilename,"Z_WAVELENGTH")
print("B coverage: {:.1f} to {:.1f} Angstroms".format(np.min(bwave),np.max(bwave)))
print("R coverage: {:.1f} to {:.1f} Angstroms".format(np.min(rwave),np.max(rwave)))
print("Z coverage: {:.1f} to {:.1f} Angstroms".format(np.min(zwave),np.max(zwave)))
Now that we understand the Data Model, let's plot some spectra. To start, let's use the file we've already been manipulating and read in the flux to go with the wavelengths we already have.
In [18]:
wave = np.hstack([bwave,rwave,zwave])
In [19]:
bflux = fitsio.read(specfilename,"B_FLUX")
rflux = fitsio.read(specfilename,"R_FLUX")
zflux = fitsio.read(specfilename,"Z_FLUX")
flux = np.hstack([bflux,rflux,zflux])
Note that the wavelength arrays are 1-D (every spectrum in the spectral file is mapped to the same binning in wavelength) but the flux array (and flux_ivar, mask etc. arrays) are 2-D, because they contain multiple spectra:
In [20]:
print(wave.shape)
print(flux.shape)
Let's plot the zeroth spectrum in this file (i.e. in this HEALPix grouping):
In [21]:
spectrum = 0
plt.plot(wave,flux[spectrum])
Out[21]:
Let's plot it again, but color-coding for the arms of the spectrograph:
In [22]:
plt.plot(bwave,bflux[spectrum],color='b')
plt.plot(rwave,rflux[spectrum],color='r')
plt.plot(zwave,zflux[spectrum],color='m')
Out[22]:
What about if we only want to plot spectra of certain target classes? The targeting information is stored in the DESI_TARGET
, BGS_TARGET
and MWS_TARGET
entries of the fibermap array:
In [23]:
fm.dtype.descr
Out[23]:
and which target corresponds to which targeting bit is stored in the desitarget mask (we imported this near the beginning of the notebook.
In [24]:
desi_mask
Out[24]:
Let's find the indexes of all standard F-stars in the spectral file:
In [25]:
stds = np.where(fm["DESI_TARGET"] & desi_mask["STD_FSTAR"])[0]
print(stds)
Where were these located on the original plate-fiber mapping?
In [26]:
plt.plot(fm["RA_TARGET"],fm["DEC_TARGET"],'b,')
plt.plot(fm["RA_TARGET"][stds],fm["DEC_TARGET"][stds],'kx')
Out[26]:
It might seem strange that there are 10 standard stars, but only 4 crosses in the plot. This is because the density of standard stars is somewhat low in this data challenge, and so some standard stars are observed multiple times. These stars will have the same TARGETID
, e.g.:
In [27]:
zipper = zip(stds,fm["TARGETID"][stds],fm["RA_TARGET"][stds],fm["DEC_TARGET"][stds])
[print("INDEX: {} TARGETID: {} RA: {} DEC: {}".format(i,tid,ra,dec)) for i,tid,ra,dec in zipper]
Out[27]:
This will not be the case for the real survey (or for subsequent iterations of the Data Challenge)!
Let's take a look at the spectra of these standard stars:
In [28]:
for i in range(len(stds)):
plt.plot(wave,flux[stds[i]],'g')
plt.show()
These seem realistic. Let's zoom in on some of the Balmer series for the zeroth standard:
In [29]:
Balmer = [4102,4341,4861,6563]
halfwindow = 50
for i in range(len(Balmer)):
plt.axis([Balmer[i]-halfwindow,Balmer[i]+halfwindow,0,np.max(flux[stds[0]])])
plt.plot(wave,flux[stds[0]])
plt.show()
The directory from which we took these spectra, also contains information on the best-fit redshifts for the spectra from the redrock code. Let's read that file and examine its contents:
In [30]:
zfilename = glob(os.path.join(basedir,'zbest*fits'))[0]
zs = fitsio.read(zfilename)
zs.dtype.descr
Out[30]:
Note that there are a different number of entries in the redshift and spectral file, meaning that there isn't a row-by-row mapping between spectra and redshifts. Ultimately, this will be fixed so that there is a row-by-row mapping between these files.
In [31]:
print(zs.shape)
print(bflux.shape)
But, the TARGETID
(which is intended to be unique) is in this file, too, allowing sources to be uniquely mapped from targeting, to spectra, to redshift. Let's extract all sources that were targeted as quasars using the fibermap information from the spectral file, and plot the first 20:
In [32]:
qsos = np.where(fm["DESI_TARGET"] & desi_mask["QSO"])[0]
for i in range(20):
plt.plot(wave,flux[qsos[i]],'b')
plt.show()
Let's match these quasar targets to the redshift file on TARGETID
to extract their best-fit redshifts from redrock
:
In [33]:
dd = defaultdict(list)
for index, item in enumerate(zs["TARGETID"]):
dd[item].append(index)
zqsos = [index for item in fm[qsos]["TARGETID"] for index in dd[item] if item in dd]
That might be hard to follow at first glance, but all I did was use some "standard" python syntax to match the indices in zs
(the ordering of objects in the redrock
redshift file) to those for quasars in fm
(the ordering of quasars in the fibermap file), on the unique TARGETID
, such that the indices stored in qsos
for fm
point to the corresponding indices in zqsos
for zs
. This might help illustrate the result:
In [34]:
zs[zqsos]["TARGETID"][0:7], fm[qsos]["TARGETID"][0:7]
Out[34]:
Let's see what best-fit template redrock
assigned to each quasar. This information is stored in the SPECTYPE
column.
In [35]:
zs[zqsos]["SPECTYPE"]
Out[35]:
Here, we've (partially) encountered a second minor bug in the Data Challenge: redrock
is currently not doing a fantastic job of fitting quasars. It correctly identified some, though. And, of course, not everything targeted as a quasar will turn out to actually be a quasar (there are some contaminants).
If we'd instead performed this same check for the standard stars, everything would have looked more reasonable:
In [36]:
dd = defaultdict(list)
for index, item in enumerate(zs["TARGETID"]):
dd[item].append(index)
zstds = [index for item in fm[stds]["TARGETID"] for index in dd[item] if item in dd]
For stars, we can also display the type of star that redrock
fit (this is stored in the SUBTYPE
column):
In [37]:
zipper = zip(zs[zstds]["SUBTYPE"],zs[zstds]["SPECTYPE"])
[ print("{}-{}".format(sub.decode('utf-8'),spec.decode('utf-8'))) for sub,spec in zipper ]
Out[37]:
(here the conversion to utf-8
is simply for display purposes because the strings in SUBTYPE
and SPECTYPE
are stored as bytes instead of unicode).
OK, back to our quasars. Let's plot the quasar targets that are identified as quasars , but add a label for the SPECTYPE
and the redshift fit by redrock
. I'll also over-plot some (approximate) typical quasar emission lines at the redrock redshift (if those lines would fall in the DESI wavelength coverage):
In [38]:
qsoid = np.where(zs[zqsos]["SPECTYPE"] == b'QSO ')[0]
#ADM bug...to be uncommented later
#qsoid = np.where(zs[zqsos]["SPECTYPE"] == b'QSO')[0]
qsolines = np.array([1216,1546,1906,2800,4853,4960,5008])
for i in range(len(qsoid)):
spectype = zs[zqsos[qsoid[i]]]["SPECTYPE"].decode('utf-8')
z = zs[zqsos[qsoid[i]]]["Z"]
plt.plot(wave,flux[qsos[qsoid[i]]],'b')
plt.suptitle("{}, z={:.3f}".format(spectype,z))
for line in qsolines:
if ((1+z)*line > np.min(wave)) & ((1+z)*line < np.max(wave)):
plt.plot([(1+z)*line,(1+z)*line],[np.min(flux[qsos[qsoid[i]]]),np.max(flux[qsos[qsoid[i]]])],'yo')
plt.show()