The goal of this notebook is to demonstrate how to generate some simple DESI spectra using the quickgen utility. For simplicity we will only generate 1D spectra and skip the more computationally intensive (yet still instructive!) step of extracting 1D spectra from simulated 2D spectra (i.e., so-called "pixel-level simulations").
For additional (albeit somewhat outdated) information and documentation about quickgen see
https://desi.lbl.gov/DocDB/cgi-bin/private/ShowDocument?docid=1429
The heart of quickgen is the SpecSim package, which you can read out here:
http://specsim.readthedocs.io/en/stable
If you identify any errors or have requests for additional functionality please create a new issue on
https://github.com/desihub/desisim/issues
or send a note to desi-data@desi.lbl.gov.
In [1]:
import os
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.table import Table
In [2]:
import desispec.io
import desisim.io
from desisim.obs import new_exposure
from desisim.scripts import quickgen
from desispec.scripts import group_spectra
In [3]:
%pylab inline
Next, make sure we have all the right environment variables set (assuming the bash shell). If any of these environment variables are missing please set them in your .bashrc file (and then restart this notebook) or create them for just this notebook session using the %set_env magic command, as we demonstrate below.
In [4]:
def check_env():
for env in ('DESIMODEL', 'DESI_ROOT', 'DESI_SPECTRO_SIM', 'DESI_SPECTRO_DATA',
'DESI_SPECTRO_REDUX', 'SPECPROD', 'PIXPROD'):
if env in os.environ:
print('{} environment set to {}'.format(env, os.getenv(env)))
else:
print('Required environment variable {} not set!'.format(env))
In [5]:
check_env()
Let's reassign the $SPECPROD environment to something other than dailytest so that we don't conflict with the outputs of the standard DESI integration test. In addition, we need to make raw data input $DESI_SPECTO_DATA match $DESI_SPECTRO_SIM/$PIXPROD where the simulated data will be written.
In [6]:
%set_env SPECPROD=example
%set_env PIXPROD=example
rawdata_dir = desisim.io.simdir()
%set_env DESI_SPECTRO_DATA=$rawdata_dir
print('Simulated raw data will be written to {}'.format(desisim.io.simdir()))
print('Pipeline will read raw data from {}'.format(desispec.io.rawdata_root()))
print(' (without knowing that it was simulated)')
print('Pipeline will write processed data to {}'.format(desispec.io.specprod_root()))
Next, let's specify the number and spectral type distribution of spectra we want to simulate, and the random seed. Setting the seed here (which can be any number at all!) ensures that your simulations are reproducible. Let's also explicitly set the night of the "observations" (the default is to use the current date) and the expid or exposure ID number (which would allow you to simulate more than one DESI exposure).
The flavor option is used to choose the correct sky-brightness model and it also determines the distribution of targets for a given flavor. For example, flavor='dark' returns the right relative sampling density of ELGs, LRGs, and QSOs. The other available (science target) options for flavor are 'dark', 'gray', 'grey', 'bright', 'bgs', 'mws', 'lrg', 'elg', 'qso', and 'std'. (You can also set flavor to either 'arc' or 'flat' but that would be boring!)
In [7]:
nspec = 100
seed = 555
flavor = 'dark'
night = '20170615'
expid = 0
The first step is to generate the fibermap and simspec files needed by quickgen. The fibermap table contains (simulated) information about the position of each target in the DESI focal plane, while the simspec table holds the "truth" spectra and the intrinsic properties of each object (redshift, noiseless photometry, [OII] flux, etc.).
In detail, the simspec and fibermap data models are described at
To generate these files we'll use new_exposure, a convenience function for generating random typical exposures of various types for testing. However, note that new_exposure isn't intended for every possible analysis; if you want to use your own mix of objects, you just need to write your own fibermap and simspec files following that format instead of calling new_exposure.
Note that in our call to new_exposure, the tileid and exptime (exposure time) optional inputs are shown for demonstration purposes but do not need to be explicitly set. In particular, the default exposure time is based on the value specified in the $DESIMODEL/data/desi.yaml parameter file.
Note: The simspec file format may change in the near future, so structure your code to separate generating spectra from the code for writing these particular formats.
In [8]:
fibermap, truth = new_exposure(flavor=flavor, nspec=nspec, seed=seed, night=night,
expid=expid, tileid=None, exptime=None)
In [9]:
rawdata_dir = desispec.io.rawdata_root()
!find $rawdata_dir | sort
Let's go a step further and read the fibermap and simspec files from on-disk.
Note that in general code should not generate filepaths by hand, but rather call desispec.io.findfile to find the file it needs. If you need to override the standard environment variable locations, you can use the outdir option, while still letting it construct the canonical filename for each type of file.
In [10]:
fiberfile = desispec.io.findfile('fibermap', night=night, expid=expid)
simspecfile = desisim.io.findfile('simspec', night=night, expid=expid)
In [11]:
print('Reading fibermap file {}'.format(fiberfile))
hdu = fits.open(fiberfile)
hdu.info()
fibermap = Table(hdu['FIBERMAP'].data)
hdu.close()
In [12]:
fibermap[:3]
Out[12]:
In [13]:
print('Reading simspec file {}.'.format(simspecfile))
hdu = fits.open(simspecfile)
hdu.info()
meta = Table(hdu['METADATA'].data)
hdu.close()
In [14]:
meta[:3]
Out[14]:
In [15]:
allobjtype = meta['OBJTYPE']
redlim = (-0.2, 1.1*meta['REDSHIFT'].max())
fig, ax = plt.subplots()
for objtype in sorted(set(allobjtype)):
indx = objtype == allobjtype
hh = ax.hist(meta['REDSHIFT'][indx], bins=nspec//3,
label=objtype, alpha=0.5, range=redlim)
ax.set_xlabel('Redshift')
ax.set_ylabel('Number of Simulated Spectra')
ax.legend(loc='upper right', ncol=3)
ax.margins(0.2)
ax.set_xlim(redlim)
Out[15]:
We're now ready to simulate DESI spectra using quickgen! Since we're calling quickgen from within this notebook (rather than from the command line in a terminal) we have to parse the simspec and fibermap filename inputs first.
quickgen generates four types of files and writes them to the \$DESI_SPECTRO_REDUX/\$SPECPROD/exposures directory: calib, sky, cframe, and frame files. We will use the cframe, or calibrated frame files, which contain the flux-calibrated and sky-subtracted DESI spectra (one file per brz camera and spectrograph).
The data model and the other files and their contents are documented here:
http://desidatamodel.readthedocs.io/en/latest/DESI_SPECTRO_REDUX/PRODNAME/exposures/NIGHT/EXPID/index.html
The code in the following cell calls the equivalent of the command line:
quickgen --simspec {simspecfile} --fibermap {fiberfile}
In [16]:
args = quickgen.parse([
'--simspec', simspecfile,
'--fibermap', fiberfile
])
quickgen.main(args)
In [17]:
cframefile = desispec.io.findfile('cframe', night=night, expid=expid, camera='b0')
print('Reading {}'.format(cframefile))
cframe = desispec.io.frame.read_frame(cframefile)
In [18]:
dir(cframe)
Out[18]:
Let's make a quick plot of the zeroth spectrum.
In [19]:
print(cframe.wave.shape, cframe.flux.shape)
In [20]:
fig, ax = plt.subplots()
ax.errorbar(cframe.wave, cframe.flux[0, :], 1/np.sqrt(cframe.ivar[0, :]))
ax.set_xlabel('Wavelength (A)')
ax.set_ylabel('Flux ($10^{-17}$ erg/s/cm$^2$)')
Out[20]:
As you can imagine, working with cframe files is pretty tedious, especially across three cameras, 10 spectrographs, and more than 35 million targets! Therefore, let's combine and reorganize the individual cframe files into spectra files grouped on the sky. Spectra are organized into healpix pixels (here chosen to have nside=64). If you're interested, you can read more about the healpix directory structure here:
https://github.com/desihub/desispec/blob/master/doc/nb/Intro_to_DESI_spectra.ipynb
Regrouping is especially important for real observations with overlapping tiles where the same object could be reobserved on different exposures separated by short or large amounts of time.
TODO: Add the spectra files to desidatamodel.
To regroup the spectra, we will run the (notebook) equivalent of the command:
desi_group_spectra --hpxnside 64
In [21]:
nside = 64
args = group_spectra.parse(['--hpxnside', '{}'.format(nside)])
group_spectra.main(args)
In [22]:
reduxdir = desispec.io.specprod_root()
!find $reduxdir | sort
As a quick example, let's plot up the zeroth spectrum in healpix pixel 19435.
In [23]:
specfilename = desispec.io.findfile('spectra', groupname=19435, nside=nside)
In [24]:
print('Reading {}'.format(specfilename))
specobj = desispec.io.read_spectra(specfilename)
In [25]:
dir(specobj)
Out[25]:
In [26]:
specobj.wave.keys(), specobj.flux.keys()
Out[26]:
In [27]:
thisone = 0
fig, ax = plt.subplots()
for camera, color in zip( ('b', 'r', 'z'), ('blue', 'red', 'magenta') ):
ax.plot(specobj.wave[camera], specobj.flux[camera][thisone], color=color)
ax.set_xlabel('Wavelength (A)')
ax.set_ylabel('Flux ($10^{-17}$ erg/s/cm$^2$)')
Out[27]: