"""
desisim.obs
===========
Utility functions related to simulating observations for DESI
"""
from __future__ import absolute_import, division, print_function
import os, sys
import numpy as np
import yaml
import sqlite3
import fcntl
import time
from astropy.io import fits
from astropy import table
import astropy.units as u
import desimodel.io
import desispec.io
from desispec.interpolation import resample_flux
from desiutil.log import get_logger
log = get_logger()
from .targets import get_targets_parallel
from . import io
import desisim.simexp
from .simexp import simulate_spectra
[docs]def new_exposure(program, nspec=5000, night=None, expid=None, tileid=None,
nproc=None, seed=None, obsconditions=None,
specify_targets=dict(), testslit=False, exptime=None,
arc_lines_filename=None, flat_spectrum_filename=None,
outdir=None, overwrite=False):
"""
Create a new exposure and output input simulation files.
Does not generate pixel-level simulations or noisy spectra.
Args:
program (str): 'arc', 'flat', 'bright', 'dark', 'bgs', 'mws', ...
nspec (int, optional): number of spectra to simulate
night (str, optional): YEARMMDD string
expid (int, optional): positive integer exposure ID
tileid (int, optional): integer tile ID
nproc (object, optional): What does this do?
seed (int, optional): random seed
obsconditions (str or dict-like, optional): see options below
specify_targets (dict of dicts, optional): Define target properties like magnitude and redshift
for each target class. Each objtype has its own key,value pair
see simspec.templates.specify_galparams_dict()
or simsepc.templates.specify_starparams_dict()
testslit (bool, optional): simulate test slit if True, default False; only for arc/flat
exptime (float, optional): exposure time [seconds], overrides obsconditions['EXPTIME']
arc_lines_filename (str, optional): use alternate arc lines filename (used if program="arc")
flat_spectrum_filename (str, optional): use alternate flat spectrum filename (used if program="flat")
outdir (str, optional): output directory
overwrite (bool, optional): optionally clobber existing files
Returns:
science: sim, fibermap, meta, obsconditions, objmeta
Writes to outdir or $DESI_SPECTRO_SIM/$PIXPROD/{night}/
* fibermap-{expid}.fits
* simspec-{expid}.fits
input obsconditions can be a string 'dark', 'gray', 'bright', or dict-like
observation metadata with keys SEEING (arcsec), EXPTIME (sec), AIRMASS,
MOONFRAC (0-1), MOONALT (deg), MOONSEP (deg). Output obsconditions is
is expanded dict-like structure.
program is used to pick the sky brightness, and is propagated to
desisim.targets.sample_objtype() to get the correct distribution of
targets for a given program, e.g. ELGs, LRGs, QSOs for program='dark'.
if program is 'arc' or 'flat', then `sim` is truth table with keys
FLUX and WAVE; and meta=None and obsconditions=None.
Also see simexp.simarc(), .simflat(), and .simscience(), the last of
which simulates a science exposure given surveysim obsconditions input,
fiber assignments, and pre-generated mock target spectra.
"""
if expid is None:
expid = get_next_expid()
if tileid is None:
tileid = get_next_tileid()
if night is None:
#- simulation obs time = now, even if sun is up
dateobs = time.gmtime()
night = get_night(utc=dateobs)
else:
#- 10pm on night YEARMMDD
night = str(night) #- just in case we got an integer instead of string
dateobs = time.strptime(night+':22', '%Y%m%d:%H')
outsimspec = desisim.io.findfile('simspec', night, expid)
outfibermap = desisim.io.findfile('simfibermap', night, expid)
if outdir is not None:
outsimspec = os.path.join(outdir, os.path.basename(outsimspec))
outfibermap = os.path.join(outdir, os.path.basename(outfibermap))
program = program.lower()
log.debug('Generating {} targets'.format(nspec))
header = dict(NIGHT=night, EXPID=expid, PROGRAM=program)
if program in ('arc', 'flat'):
header['FLAVOR'] = program
else:
header['FLAVOR'] = 'science'
#- ISO 8601 DATE-OBS year-mm-ddThh:mm:ss
header['DATE-OBS'] = time.strftime('%FT%T', dateobs)
if program == 'arc':
if arc_lines_filename is None :
infile = os.getenv('DESI_ROOT')+'/spectro/templates/calib/v0.4/arc-lines-average-in-vacuum-from-winlight-20170118.fits'
else :
infile = arc_lines_filename
arcdata = fits.getdata(infile, 1)
#- clip unphysical negative values in arc template
arcdata['ELECTRONS'] = arcdata['ELECTRONS'].clip(0)
if exptime is None:
exptime = 5
wave, phot, fibermap = desisim.simexp.simarc(arcdata, nspec=nspec, testslit=testslit)
header['EXPTIME'] = exptime
desisim.io.write_simspec_arc(outsimspec, wave, phot, header, fibermap=fibermap, overwrite=overwrite)
fibermap.meta['NIGHT'] = night
fibermap.meta['EXPID'] = expid
desispec.io.write_fibermap(outfibermap, fibermap)
truth = dict(WAVE=wave, PHOT=phot, UNITS='photon')
return truth, fibermap, None, None, None
elif program == 'flat':
if flat_spectrum_filename is None :
infile = os.getenv('DESI_ROOT')+'/spectro/templates/calib/v0.4/flat-3100K-quartz-iodine.fits'
else :
infile = flat_spectrum_filename
if exptime is None:
exptime = 10
sim, fibermap = desisim.simexp.simflat(infile, nspec=nspec,
exptime=exptime, testslit=testslit, psfconvolve=False)
header['EXPTIME'] = exptime
header['FLAVOR'] = 'flat'
desisim.io.write_simspec(sim, truth=None, fibermap=fibermap, obs=None,
expid=expid, night=night, header=header, filename=outsimspec, overwrite=overwrite)
fibermap.meta['NIGHT'] = night
fibermap.meta['EXPID'] = expid
desispec.io.write_fibermap(outfibermap, fibermap)
# fluxunits = 1e-17 * u.erg / (u.s * u.cm**2 * u.Angstrom)
fluxunits = '1e-17 erg/(s * cm2 * Angstrom)'
flux = sim.simulated['source_flux'].to(fluxunits)
wave = sim.simulated['wavelength'].to('Angstrom')
truth = dict(WAVE=wave, FLUX=flux, UNITS=str(fluxunits))
return truth, fibermap, None, None, None
#- all other programs
fibermap, (flux, wave, meta, objmeta) = get_targets_parallel(nspec, program,
tileid=tileid, nproc=nproc, seed=seed, specify_targets=specify_targets)
fibermap["PETAL_LOC"] = fibermap["FIBER"]//500
if obsconditions is None:
if program in ['dark', 'lrg', 'qso']:
obsconditions = desisim.simexp.reference_conditions['DARK']
elif program in ['elg', 'gray', 'grey']:
obsconditions = desisim.simexp.reference_conditions['GRAY']
elif program in ['mws', 'bgs', 'bright']:
obsconditions = desisim.simexp.reference_conditions['BRIGHT']
else:
raise ValueError('unknown program {}'.format(program))
elif isinstance(obsconditions, str):
try:
obsconditions = desisim.simexp.reference_conditions[obsconditions.upper()]
except KeyError:
raise ValueError('obsconditions {} not in {}'.format(
obsconditions.upper(),
list(desisim.simexp.reference_conditions.keys())))
if exptime is not None:
obsconditions['EXPTIME'] = exptime
sim = simulate_spectra(wave, flux, fibermap=fibermap,
obsconditions=obsconditions, psfconvolve=False)
#- Write fibermap
telera, teledec = io.get_tile_radec(tileid)
hdr = dict(
NIGHT = (night, 'Night of observation YEARMMDD'),
EXPID = (expid, 'DESI exposure ID'),
TILEID = (tileid, 'DESI tile ID'),
PROGRAM = (program, 'program [dark, bright, ...]'),
FLAVOR = ('science', 'Flavor [arc, flat, science, zero, ...]'),
TELRA = (telera, 'Telescope pointing RA [degrees]'),
TELDEC = (teledec, 'Telescope pointing dec [degrees]'),
AIRMASS = (obsconditions['AIRMASS'], 'Airmass at middle of exposure'),
EXPTIME = (obsconditions['EXPTIME'], 'Exposure time [sec]'),
SEEING = (obsconditions['SEEING'], 'Seeing FWHM [arcsec]'),
MOONFRAC = (obsconditions['MOONFRAC'], 'Moon illumination fraction 0-1; 1=full'),
MOONALT = (obsconditions['MOONALT'], 'Moon altitude [degrees]'),
MOONSEP = (obsconditions['MOONSEP'], 'Moon:tile separation angle [degrees]'),
)
hdr['DATE-OBS'] = (time.strftime('%FT%T', dateobs), 'Start of exposure')
simfile = io.write_simspec(sim, meta, fibermap, obsconditions,
expid, night, objmeta=objmeta, header=hdr,
filename=outsimspec, overwrite=overwrite)
if not isinstance(fibermap, table.Table):
fibermap = table.Table(fibermap)
fibermap.meta.update(hdr)
desispec.io.write_fibermap(outfibermap, fibermap)
log.info('Wrote '+outfibermap)
update_obslog(obstype='science', program=program, expid=expid, dateobs=dateobs, tileid=tileid)
return sim, fibermap, meta, obsconditions, objmeta
#- Mapping of DESI objtype to things specter knows about
[docs]def specter_objtype(desitype):
'''
Convert a list of DESI object types into ones that specter knows about
'''
intype = np.atleast_1d(desitype)
desi2specter = dict(
STAR='STAR', STD='STAR', MWS_STAR='STAR',
LRG='LRG', ELG='ELG', QSO='QSO', QSO_BAD='STAR',
BGS='LRG', # !!!
SKY='SKY'
)
unknown_types = set(intype) - set(desi2specter.keys())
if len(unknown_types) > 0:
raise ValueError('Unknown input objtypes {}'.format(unknown_types))
results = np.zeros(len(intype), dtype=(str, 8))
for objtype in desi2specter:
ii = (intype == objtype)
results[ii] = desi2specter[objtype]
assert np.count_nonzero(results == '') == 0
if isinstance(desitype, str):
return results[0]
else:
return results
[docs]def get_next_tileid(program='DARK'):
"""
Return tileid of next tile to observe
Args:
program (optional): dark, gray, or bright
Note:
Simultaneous calls will return the same tileid;
it does *not* reserve the tileid.
"""
program = program.upper()
if program not in ('DARK', 'GRAY', 'GREY', 'BRIGHT',
'ELG', 'LRG', 'QSO', 'LYA', 'BGS', 'MWS'):
return -1
#- Read DESI tiling and trim to just tiles in DESI footprint
tiles = table.Table(desimodel.io.load_tiles())
#- HACK: update tilelist to include PROGRAM, etc.
if 'PROGRAM' not in tiles.colnames:
log.error('You are using an out-of-date desi-tiles.fits file from desimodel')
log.error('please update your copy of desimodel/data')
log.warning('proceeding anyway with a workaround for now...')
tiles['PASS'] -= min(tiles['PASS']) #- standardize to starting at 0 not 1
brighttiles = tiles[tiles['PASS'] <= 2].copy()
brighttiles['TILEID'] += 50000
brighttiles['PASS'] += 5
tiles = table.vstack([tiles, brighttiles])
program_col = table.Column(name='PROGRAM', length=len(tiles), dtype=(str, 6))
tiles.add_column(program_col)
tiles['PROGRAM'][tiles['PASS'] <= 3] = 'DARK'
tiles['PROGRAM'][tiles['PASS'] == 4] = 'GRAY'
tiles['PROGRAM'][tiles['PASS'] >= 5] = 'BRIGHT'
else:
tiles['PROGRAM'] = np.char.strip(tiles['PROGRAM'])
#- If obslog doesn't exist yet, start at tile 0
dbfile = io.simdir()+'/etc/obslog.sqlite'
if not os.path.exists(dbfile):
obstiles = set()
else:
#- Read obslog to get tiles that have already been observed
db = sqlite3.connect(dbfile)
result = db.execute('SELECT tileid FROM obslog WHERE program="{}"'.format(program))
obstiles = set( [row[0] for row in result] )
db.close()
#- Just pick the next tile in sequential order
program_tiles = tiles['TILEID'][tiles['PROGRAM'] == program]
nexttile = int(min(set(program_tiles) - obstiles))
log.debug('{} tiles in program {}'.format(len(program_tiles), program))
log.debug('{} observed tiles'.format(len(obstiles)))
return nexttile
[docs]def get_next_expid(n=None):
"""
Return the next exposure ID to use from {proddir}/etc/next_expid.txt
and update the exposure ID in that file.
Use file locking to prevent multiple readers from getting the same
ID or accidentally clobbering each other while writing.
Args:
n (int, optional): number of contiguous expids to return as a list.
If None, return a scalar. Note that n=1 returns a list of length 1.
BUGS:
* if etc/next_expid.txt doesn't exist, initial file creation is
probably not threadsafe.
* File locking mechanism doesn't work on NERSC Edison, to turned off
for now.
"""
#- Full path to next_expid.txt file
filename = io.simdir()+'/etc/next_expid.txt'
if not os.path.exists(io.simdir()+'/etc/'):
os.makedirs(io.simdir()+'/etc/')
#- Create file if needed; is this threadsafe? Probably not.
if not os.path.exists(filename):
fw = open(filename, 'w')
fw.write("0\n")
fw.close()
#- Open the file, but get exclusive lock before reading
f0 = open(filename)
### fcntl.flock(f0, fcntl.LOCK_EX)
expid = int(f0.readline())
#- Write update expid to the file
fw = open(filename, 'w')
if n is None:
fw.write(str(expid+1)+'\n')
else:
fw.write(str(expid+n)+'\n')
fw.close()
#- Release the file lock
### fcntl.flock(f0, fcntl.LOCK_UN)
f0.close()
if n is None:
return expid
else:
return list(range(expid, expid+n))
[docs]def get_night(t=None, utc=None):
"""
Return YEARMMDD for tonight. The night roles over at local noon.
i.e. 1am and 11am is the previous date; 1pm is the current date.
Args:
t : local time.struct_time tuple of integers
(year, month, day, hour, min, sec, weekday, dayofyear, DSTflag)
default is time.localtime(), i.e. now
utc : time.struct_time tuple for UTC instead of localtime
Note:
this only has one second accuracy; good enough for sims but *not*
to be used for actual DESI ops.
"""
#- convert t to localtime or fetch localtime if needed
if utc is not None:
t = time.localtime(time.mktime(utc) - time.timezone)
elif t is None:
t = time.localtime()
#- what date/time was it 12 hours ago? "Night" rolls over at noon local
night = time.localtime(time.mktime(t) - 12*3600)
#- format that as YEARMMDD
return time.strftime('%Y%m%d', night)
#- I'm not really sure this is a good idea.
#- I'm sure I will want to change the schema later...
[docs]def update_obslog(obstype='science', program='DARK', expid=None, dateobs=None,
tileid=-1, ra=None, dec=None):
"""
Update obslog with a new exposure
obstype : 'arc', 'flat', 'bias', 'test', 'science', ...
program : 'DARK', 'GRAY', 'BRIGHT', 'CALIB'
expid : integer exposure ID, default from get_next_expid()
dateobs : time.struct_time tuple; default time.localtime()
tileid : integer TileID, default -1, i.e. not a DESI tile
ra, dec : float (ra, dec) coordinates, default tile ra,dec or (0,0)
returns tuple (expid, dateobs)
TODO: normalize obstype vs. program; see desisim issue #97
"""
#- Connect to sqlite database file and create DB if needed
dbdir = io.simdir() + '/etc'
if not os.path.exists(dbdir):
os.makedirs(dbdir)
dbfile = dbdir+'/obslog.sqlite'
with sqlite3.connect(dbfile, isolation_level="EXCLUSIVE") as db:
db.execute("""\
CREATE TABLE IF NOT EXISTS obslog (
expid INTEGER PRIMARY KEY,
dateobs DATETIME, -- seconds since Unix Epoch (1970)
night TEXT, -- YEARMMDD
obstype TEXT DEFAULT "science",
program TEXT DEFAULT "DARK",
tileid INTEGER DEFAULT -1,
ra REAL DEFAULT 0.0,
dec REAL DEFAULT 0.0
)
""")
#- Fill in defaults
if expid is None:
expid = get_next_expid()
if dateobs is None:
dateobs = time.localtime()
if ra is None:
assert (dec is None)
if tileid < 0:
ra, dec = (0.0, 0.0)
else:
ra, dec = io.get_tile_radec(tileid)
night = get_night(utc=dateobs)
insert = """\
INSERT OR REPLACE INTO obslog(expid,dateobs,night,obstype,program,tileid,ra,dec)
VALUES (?,?,?,?,?,?,?,?)
"""
db.execute(insert, (int(expid), time.mktime(dateobs), str(night), str(obstype.upper()), str(program.upper()), int(tileid), float(ra), float(dec)))
db.commit()
return expid, dateobs