Add function to add calibration exposures

doc/changes.rst

@@ -5,7 +5,8 @@ surveysim change log
 0.8.3 (unreleased)
 ------------------
 
-* No changes yet.
+* Add ``surveysim.util.add_calibration_exposures()``, to add simulated
+  calibration exposures to a set of science exposures.
 
 0.8.2 (2017-10-09)
 ------------------

py/surveysim/test/test_util.py

@@ -0,0 +1,64 @@
+"""Test surveysim.util.
+"""
+import unittest
+import numpy as np
+from astropy.table import Table, Column
+from ..util import add_calibration_exposures
+
+class TestUtil(unittest.TestCase):
+    """Test surveysim.util.
+    """
+
+    @classmethod
+    def setUpClass(cls):
+        pass
+
+    @classmethod
+    def tearDownClass(cls):
+        pass
+
+    def setUp(self):
+        pass
+
+    def tearDown(self):
+        pass
+
+    def test_add_calibration_exposures(self):
+        """Test adding calibration exposures to science exposures.
+        """
+        exposures = Table()
+        exposures['TILEID'] = Column(np.array([0, 1], dtype=np.int32))
+        exposures['PASS'] = Column(np.array([0, 0], dtype=np.int16))
+        exposures['RA'] = Column(np.array([0.0, 1.0], dtype=np.float64))
+        exposures['DEC'] = Column(np.array([0.0, 1.0], dtype=np.float64))
+        exposures['EBMV'] = Column(np.array([0.0, 1.0], dtype=np.float64))
+        exposures['NIGHT'] = Column(np.array(['20200101', '20200102'], dtype=(str, 8)))
+        exposures['MJD'] = Column(np.array([58849.0, 58850.0], dtype=np.float64))
+        exposures['EXPTIME'] = Column(np.array([0.0, 1.0], dtype=np.float64), unit='s')
+        exposures['SEEING'] = Column(np.array([0.0, 1.0], dtype=np.float64), unit='arcsec')
+        exposures['TRANSPARENCY'] = Column(np.array([0.0, 1.0], dtype=np.float64))
+        exposures['AIRMASS'] = Column(np.array([0.0, 1.0], dtype=np.float64))
+        exposures['MOONFRAC'] = Column(np.array([0.0, 1.0], dtype=np.float64))
+        exposures['MOONALT'] = Column(np.array([0.0, 1.0], dtype=np.float64), unit='deg')
+        exposures['MOONSEP'] = Column(np.array([0.0, 1.0], dtype=np.float64), unit='deg')
+        exposures['PROGRAM'] = Column(np.array(['DARK', 'DARK'], dtype=(str, 6)))
+        exposures['FLAVOR'] = Column(np.array(['science', 'science'], dtype=(str, 7)))
+        output = add_calibration_exposures(exposures)
+        self.assertEqual(len(output), 14)
+        self.assertEqual('EXPID', output.colnames[0])
+        self.assertTrue(np.all(output['EXPID'] == np.arange(14, dtype=np.int32)))
+        self.assertEqual(output['RA'].unit, 'deg')
+        self.assertTrue(np.all(np.diff(output['MJD']) >= 0))
+        self.assertTrue(np.all(np.diff(output['EXPID']) == 1))
+        bad_exposures = Table()
+        bad_exposures['MJD'] = Column(np.array([58851.0, 58850.0], dtype=np.float64))
+        with self.assertRaises(ValueError):
+            output = add_calibration_exposures(bad_exposures)
+
+
+def test_suite():
+    """Allows testing of only this module with the command::
+
+        python setup.py test -m <modulename>
+    """
+    return unittest.defaultTestLoader.loadTestsFromName(__name__)

py/surveysim/util.py

@@ -0,0 +1,100 @@
+"""Simulation utilities that may be used by other packages.
+"""
+from __future__ import print_function, division, absolute_import
+import numpy as np
+from astropy.table import Column
+
+
+def add_calibration_exposures(exposures, flats_per_night=3, arcs_per_night=3,
+                              darks_per_night=0, zeroes_per_night=0,
+                              exptime=None, readout=30.0):
+    """Adds calibration exposures to a set of science exposures and
+    assigns exposure IDs.
+
+    Parameters
+    ----------
+    exposures : :class:`astropy.table.Table`
+        A table of science exposures, produced by *e.g.*
+        ``desisurvey.progress.Progress.get_exposures()``.
+        The exposures must be sorted by increasing MJD/timestamp.
+    flats_per_night : :class:`int`, optional
+        Add this many arc exposures per night (default 3).
+    arcs_per_night : :class:`int`, optional
+        Add this many arc exposures per night (default 3).
+    darks_per_night : :class:`int`, optional
+        Add this many dark exposures per night (default 0).
+    zeroes_per_night : :class:`int`, optional
+        Add this many zero exposures per night (default 0).
+    exptime : :class:`dict`, optional
+        A dictionary setting calibration exposure times for each
+        calibration flavor.
+    readout : :class:`float`, optional
+        Set readout time for calibration exposures (default 30.0 s).
+
+    Returns
+    -------
+    :class:`astropy.table.Table`
+        A table augmented with calibration exposures.
+
+    Raises
+    ------
+    ValueError
+        If the input is not sorted by increasing MJD/timestamp.
+    """
+    if not np.all(np.diff(exposures['MJD']) >= 0):
+        raise ValueError("Input is not sorted by increasing MJD!")
+    if exptime is None:
+        exptime = {'flat': 10.0, 'arc': 10.0, 'dark': 1000.0, 'zero': 0.0}
+    output = exposures[:0].copy()
+    expid_in_exposures = 'EXPID' in exposures.colnames
+    for c in ('RA', 'DEC'):
+        if output[c].unit is None:
+            output[c].unit = 'deg'
+    current_night = '19730703'
+    expid = 0
+    calib_time = lambda x: exptime[x] + readout
+    calib_sequence = (['flat']*flats_per_night + ['arc']*arcs_per_night +
+                      ['dark']*darks_per_night + ['zero']*zeroes_per_night)
+    calib_times = np.cumsum(np.array([calib_time(c)
+                                      for c in calib_sequence]))[::-1]
+    # None fields get filled in for each calibration exposure.
+    calib_data = {'TILEID': -1,
+                  'PASS': -1,
+                  'RA': 0.0,
+                  'DEC': 0.0,
+                  'EBMV': 0.0,
+                  'NIGHT': None,
+                  'MJD': None,
+                  'EXPTIME': None,
+                  'SEEING': 0.0,
+                  'TRANSPARENCY': 0.0,
+                  'AIRMASS': 0.0,
+                  'MOONFRAC': 0.0,
+                  'MOONALT': 0.0,
+                  'MOONSEP': 0.0,
+                  'PROGRAM': 'CALIB',
+                  'FLAVOR': None
+                  }
+    if expid_in_exposures:
+        calib_data['EXPID'] = None
+    for i, night in enumerate(exposures['NIGHT']):
+        if night != current_night:
+            current_night = night
+            for j, c in enumerate(calib_sequence):
+                if expid_in_exposures:
+                    calib_data['EXPID'] = expid
+                calib_data['NIGHT'] = night
+                calib_data['MJD'] = exposures['MJD'][i] - calib_times[j]/86400.0
+                calib_data['EXPTIME'] = exptime[c]
+                calib_data['FLAVOR'] = c
+                output.add_row(calib_data)
+                expid += 1
+        output.add_row(exposures[i])
+        if expid_in_exposures:
+            output['EXPID'][expid] = expid
+        expid += 1
+    if not expid_in_exposures:
+        output.add_column(Column(np.arange(len(output), dtype=np.int32),
+                                 name='EXPID', description='Exposure ID'),
+                          index=0)
+    return output