In [57]:
import os
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.table import Table
import yaml
In [58]:
import desispec.io
import desisim.io
from desisim.scripts import quickgen
from desispec.scripts import group_spectra
from desispec.io.util import write_bintable, makepath
In [59]:
#%pylab inline
from desiutil.log import get_logger
log = get_logger(level='WARNING')
In [60]:
# currently in the bgssim branch of desisim
from desisim.obs import new_exposure
# new_exposure(flavor, nspec=5000, night=None, expid=None, tileid=None,
# airmass=1.0, exptime=None, seed=None, testslit=False,
# arc_lines_filename=None, flat_spectrum_filename=None,
# target_densities = {})
# returns : fibermap, truth
# or
# wavelengths = qsim.source.wavelength_out.to(u.Angstrom).value
# bgs = desisim.templates.BGS(wave=wavelengths, add_SNeIa=args.add_SNeIa)
# flux, tmpwave, meta1 = bgs.make_templates(nmodel=nobj, seed=args.seed, zrange=args.zrange_bgs,
# rmagrange=args.rmagrange_bgs,sne_rfluxratiorange=args.sne_rfluxratiorange)
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 [61]:
nspec = 200
seed = 555
night = '20170701'
flavor = 'bgs'
nexp = 10 # number of exposures
In [62]:
exptime_range = (300, 300)
airmass_range = (1.25, 1.25)
moonphase_range = (0.0, 1.0)
moonangle_range = (0, 150)
moonzenith_range = (0, 60)
In [63]:
targ_dens = {}
targ_dens['frac_std'] = 0.02
targ_dens['frac_sky'] = 0.08
targ_dens['area'] = 14000.0
targ_dens['area_bgs'] = 14000
targ_dens['nobs_bgs_bright'] = 762
targ_dens['nobs_bgs_faint'] = 475
targ_dens['ntarget_bgs_bright'] = 818
targ_dens['ntarget_bgs_faint'] = 618
targ_dens['success_bgs_bright'] = 0.97
targ_dens['success_bgs_faint'] = 0.92
targ_dens['nobs_mws'] = 700
targ_dens['ntarget_mws'] = 736
targ_dens['success_mws'] = 0.99
In [64]:
def check_env():
for env in ('DESIMODEL', 'DESI_ROOT', 'DESI_SPECTRO_SIM', 'DESI_SPECTRO_DATA',
'DESI_SPECTRO_REDUX', 'SPECPROD', 'PIXPROD','DESI_BASIS_TEMPLATES'):
if env in os.environ:
print('{} environment set to {}'.format(env, os.getenv(env)))
else:
print('Required environment variable {} not set!'.format(env))
In [65]:
check_env()
In [66]:
%set_env SPECPROD=bgs-specsim-paper-kr
%set_env PIXPROD=bgs-specsim-paper-kr
rawdata_dir = desisim.io.simdir()
%set_env DESI_SPECTRO_DATA=$rawdata_dir
%set_env DESI_BASIS_TEMPLATES=/project/projectdirs/desi/spectro/templates/basis_templates/trunk
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()))
In [67]:
rand = np.random.RandomState(seed)
In [68]:
expids = np.arange(nexp).astype(int)
exptime = rand.uniform(exptime_range[0], exptime_range[1], nexp)
airmass = rand.uniform(airmass_range[0], airmass_range[1], nexp)
moonphase = rand.uniform(moonphase_range[0], moonphase_range[1], nexp)
moonangle = rand.uniform(moonangle_range[0], moonangle_range[1], nexp)
moonzenith = rand.uniform(moonzenith_range[0], moonzenith_range[1], nexp)
In [69]:
metafile = os.path.join( desisim.io.simdir(), 'mysim.fits')
metacols = [
('BRICKNAME', 'S20'),
('SEED', 'S20'),
('EXPTIME', 'f4'),
('AIRMASS', 'f4'),
('MOONPHASE', 'f4'),
('MOONANGLE', 'f4'),
('MOONZENITH', 'f4')]
meta = Table(np.zeros(nexp, dtype=metacols))
meta['EXPTIME'].unit = 's'
meta['MOONANGLE'].unit = 'deg'
meta['MOONZENITH'].unit = 'deg'
#meta['BRICKNAME'] = ['{}-{:03d}'.format(args.brickname, ii) for ii in range(args.nbrick)]
meta['EXPTIME'] = exptime
meta['AIRMASS'] = airmass
meta['MOONPHASE'] = moonphase
meta['MOONANGLE'] = moonangle
meta['MOONZENITH'] = moonzenith
log.info('Writing {}'.format(metafile))
In [70]:
targetyaml = os.path.join(os.environ['DESIMODEL'],'data','targets','targets.yaml')
tgt = yaml.load(open(targetyaml))
for key, val in targ_dens.items():
tgt[key] = val
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
In [ ]:
for ii,expid in enumerate(expids):
fibermap, truth = new_exposure(flavor, nspec=nspec, night=night, expid=int(expid), tileid=None,\
airmass=airmass[ii], exptime=exptime[ii], seed=seed,\
target_densities=tgt)
quickgenTo get around the fact that we aren't using the command line, we use the arg parser and pass the arguments to the main function of quickgen directly.
more information at: http://desidatamodel.readthedocs.io/en/latest/DESI_SPECTRO_REDUX/PRODNAME/exposures/NIGHT/EXPID/index.html
Quickgen additional commands for 'quickbrick mode:'
'--objtype', 'BGS',
'--brickname', 'whatever',
'--zrange-bgs', (0.01, 0.4),
'--rmagrange-bgs', (15.0,19.5)
'--exptime', None
'--airmass', 1.5
In [ ]:
for ii,expid in enumerate(expids):
fiberfile = desispec.io.findfile('fibermap', night=night, expid=expid)
simspecfile = desisim.io.findfile('simspec', night=night, expid=expid)
args = quickgen.parse([
'--simspec', simspecfile,
'--fibermap', fiberfile,
'--nspec', str(nspec),
'--seed', str(seed),
'--moon-phase', str(moonphase[ii]),
'--moon-angle', str(moonangle[ii]),
'--moon-zenith', str(moonzenith[ii])
])
quickgen.main(args)
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.
In [86]:
nside = 64
args = group_spectra.parse(['--hpxnside', '{}'.format(nside)])
group_spectra.main(args)