Merge pull request #473 from LSSTDESC/u/jchiang/header_fixes

Header fixes and LSSTCamSim support

.gitattributes

@@ -0,0 +1,3 @@
+*.fits binary
+*.fits.fz binary
+*.fits.gz binary

data/LsstCamSim_info.yaml

@@ -0,0 +1,6 @@
+tree_rings_file_name: tree_ring_parameters_2018-04-26.txt
+vignetting_file_name: LSSTCam_vignetting_data.json
+telescope_format: LSST_%s.yaml
+bias_levels_file: LSSTCam_bias_levels_run_13421.json
+camera_name: LsstCam
+ndets: 189

imsim/bandpass.py

@@ -89,7 +89,7 @@ def RubinBandpass(
     if (camera is None) != (det_name is None):
         raise ValueError("Must provide both camera and det_name if using one.")
     match camera:
-        case 'LsstCam':
+        case 'LsstCam' | 'LsstCamSim' :
             camera = 'lsstCam'
         case 'LsstComCamSim':
             camera = 'comCamSim'

imsim/camera.py

@@ -162,7 +162,7 @@ def get_camera(camera='LsstCam'):
     -------
     lsst.afw.cameraGeom.Camera
     """
-    valid_cameras = ('LsstCam', 'LsstCamImSim', 'LsstComCamSim')
+    valid_cameras = ('LsstCam', 'LsstCamSim', 'LsstCamImSim', 'LsstComCamSim')
     if camera not in valid_cameras:
         raise ValueError('Invalid camera: %s', camera)
     if camera not in _camera_cache:

imsim/ccd.py

@@ -170,9 +170,11 @@ def parse(item, type, default):
         # Now construct the image header
         image.header['MJD'] = mjd
         image.header['MJD-OBS'] = mjd_obs, 'Start of exposure'
-        # NOTE: Should this day be the current day,
-        # or the day at the time of the most recent noon?
-        dayobs = astropy.time.Time(mjd_obs, format='mjd').strftime('%Y%m%d')
+        # Subtract half-day offset from mjd-obs for dayobs calculation
+        # following Rubin convention.  See
+        # https://github.com/lsst/astro_metadata_translator/blob/w.2024.19/python/astro_metadata_translator/translator.py#L1065
+        # and Appendix A of https://docushare.lsst.org/docushare/dsweb/Get/LSE-400
+        dayobs = astropy.time.Time(mjd_obs - 0.5, format='mjd').strftime('%Y%m%d')
         image.header['DAYOBS'] = dayobs
         image.header['SEQNUM'] = seqnum
         image.header['CONTRLLR'] = 'S', 'simulated data'
@@ -192,6 +194,8 @@ def parse(item, type, default):
         image.header['AMSTART'] = airmass
         image.header['AMEND'] = airmass  # wrong, does anyone care?
         image.header['FOCUSZ'] = parse('focusZ', float, 0.0)
+        image.header['ALTITUDE'] = parse('altitude', float, 'N/A')
+        image.header['AZIMUTH'] = parse('azimuth', float, 'N/A')
 
         # If there's anything left in header_vals, add it to the header.
         for k in header_vals:

imsim/opsim_data.py

@@ -114,8 +114,9 @@ def _read_opsim_db(self):
                 self.meta[key] = value
 
             # Determine the daily sequence number for this exposure by
-            # counting the number of snaps since int(observationStartMJD).
-            t0 = int(self.meta['observationStartMJD'])
+            # counting the number of snaps since int(observationStartMJD)
+            # with a half-day offset to account for Rubin's dayobs definition.
+            t0 = int(self.meta['observationStartMJD']) - 0.5
             sql = f'''select numExposures from observations where
                       {t0} <= observationStartMJD and
                       observationId < {self.visit}'''

imsim/readout.py

@@ -259,7 +259,6 @@ def get_primary_hdu(eimage, lsst_num, camera_name=None,
         telcode = 'MC'
     elif camera_name == 'LsstComCamSim' :
         phdu.header['FILTER'] = ComCam_filter_map.get(band, None)
-        phdu.header['TELESCOP'] = SIMONYI_TELESCOPE
         phdu.header['INSTRUME'] = 'ComCamSim'
         phdu.header['RAFTBAY'] = raft
         phdu.header['CCDSLOT'] = sensor
@@ -270,16 +269,31 @@ def get_primary_hdu(eimage, lsst_num, camera_name=None,
         telcode = 'CC'
     else:
         phdu.header['FILTER'] = LSSTCam_filter_map.get(band, None)
-        phdu.header['INSTRUME'] = 'LSSTCam'
+        if camera_name == 'LsstCamSim':
+            phdu.header['INSTRUME'] = 'LSSTCamSim'
+        else:
+            phdu.header['INSTRUME'] = 'LSSTCam'
         phdu.header['RAFTBAY'] = raft
         phdu.header['CCDSLOT'] = sensor
         phdu.header['RA'] = ratel
         phdu.header['DEC'] = dectel
         phdu.header['ROTCOORD'] = 'sky'
+        phdu.header['ROTPA'] = rotang
         telcode = 'MC'
     dayobs = eimage.header['DAYOBS']
     seqnum = eimage.header['SEQNUM']
     contrllr = eimage.header['CONTRLLR']
+    phdu.header['TELESCOP'] = SIMONYI_TELESCOPE
+    phdu.header['TELCODE'] = telcode
+    phdu.header['RASTART'] = ratel
+    phdu.header['DECSTART'] = dectel
+    phdu.header['ELSTART'] = eimage.header['ALTITUDE']
+    phdu.header['AZSTART'] = eimage.header['AZIMUTH']
+    if eimage.header['IMGTYPE'] == 'SKYEXP':
+        phdu.header['RADESYS'] = 'ICRS'
+        phdu.header['TRACKSYS'] = 'RADEC'
+    else:
+        phdu.header['TRACKSYS'] = 'LOCAL'
     phdu.header['OBSID'] = f"{telcode}_{contrllr}_{dayobs}_{seqnum:06d}"
     phdu.header['MJD-OBS'] = mjd_obs
     phdu.header['HASTART'] = eimage.header['HASTART']
@@ -288,9 +302,10 @@ def get_primary_hdu(eimage, lsst_num, camera_name=None,
     phdu.header['DATE-END'] = Time(mjd_end, format='mjd', scale='tai').to_value('isot')
     phdu.header['AMSTART'] = eimage.header['AMSTART']
     phdu.header['AMEND'] = eimage.header['AMEND']
+    phdu.header['ORIGIN'] = "imSim"
     phdu.header['IMSIMVER'] = __version__
     phdu.header['PKG00000'] = 'throughputs'
-    phdu.header['VER00000'] = '1.4'
+    phdu.header['VER00000'] = '1.9'
     phdu.header['CHIPID'] = det_name
     phdu.header['FOCUSZ'] = eimage.header['FOCUSZ']
 

tests/test_object_positions.py

@@ -114,6 +114,9 @@ def test_output_catalog():
     assert data['realized_flux'] > 0.99 * data['phot_flux']
     # fft_flux is 0 when object was photon shot.
     assert data['fft_flux'] == 0.
+    output_dir = "fits_LsstCam"
+    if os.path.isdir(output_dir):
+        shutil.rmtree(output_dir)
 
 
 if __name__ == "__main__":

tests/test_raw_file_writing.py

@@ -29,7 +29,7 @@ def test_raw_file_writing(self):
         This is mostly an operational test that the raw files can be
         written from a galsim image.
         """
-        eimage_file = str(DATA_DIR / 'eimage_00449053-1-r-R22_S11-det094.fits')
+        eimage_file = str(DATA_DIR / 'eimage_00449053-1-r-R22_S11-det094.fits.gz')
         eimage = galsim.fits.read(eimage_file)
         eimage.header = galsim.FitsHeader(file_name=eimage_file)
 

tests/test_readout.py

@@ -20,7 +20,7 @@ class ImageSourceTestCase(unittest.TestCase):
     "TestCase class for ImageSource."
 
     def setUp(self):
-        self.eimage_file = str(DATA_DIR / 'eimage_00449053-1-r-R22_S11-det094.fits')
+        self.eimage_file = str(DATA_DIR / 'eimage_00449053-1-r-R22_S11-det094.fits.gz')
         instcat_file = str(DATA_DIR / 'tiny_instcat.txt')
         self.image = galsim.fits.read(self.eimage_file)
         self.image.header = galsim.FitsHeader(file_name=self.eimage_file)
@@ -94,11 +94,31 @@ def test_raw_file_headers(self):
         "Test contents of raw file headers."
         outfile = 'raw_file_test.fits'
         self.readout.writeFile(outfile, self.readout_config, self.config, self.logger)
+        # Some required keywords that were noted as missing for OR3.
+        # See https://github.com/LSSTDESC/imSim/issues/457.
+        expected_keywords = [
+            "RA",
+            "DEC",
+            "RASTART",
+            "DECSTART",
+            "ROTPA",
+            "ROTCOORD",
+            "HASTART",
+            "ELSTART",
+            "AZSTART",
+            "AMSTART",
+            "TRACKSYS",
+            "RADESYS",
+            "ORIGIN",
+            "TELCODE"
+        ]
         with fits.open(outfile) as hdus:
             self.assertEqual(hdus[0].header['IMSIMVER'], imsim.__version__)
             # Test added_keywords are included correctly
             self.assertEqual(hdus[0].header['TESTKEY1'], 'TESTVAL1')
             self.assertEqual(hdus[0].header['SOMEMATH'], '3')
+            for keyword in expected_keywords:
+                hdus[0].header[keyword]
         os.remove(outfile)
 
     def test_no_opsim(self):