-
Notifications
You must be signed in to change notification settings - Fork 22
/
obs.py
457 lines (374 loc) · 16.9 KB
/
obs.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
"""
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
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)
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
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
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
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))
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...
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