Source code for desisim.obs

"""
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