Merge f66a431 into 08eb943

bin.src/write_amp_image_file.py

@@ -16,14 +16,20 @@
                     default='INFO', help='Logging level. Default: "INFO"')
 parser.add_argument("--opsim_db", default=None, type=str,
                     help="OpSim db file as alternative source of pointing info")
-args = parser.parse_args()
+parser.add_argument('--config_file', type=str, default=None,
+                    help="Config file. If None, the default config will be used.")
 
-desc.imsim.read_config()
+args = parser.parse_args()
 
 logger = desc.imsim.get_logger(args.log_level)
 
+config = desc.imsim.read_config(args.config_file)
+
 image_source = desc.imsim.ImageSource.create_from_eimage(args.eimage_file,
                                                          opsim_db=args.opsim_db,
                                                          logger=logger)
-outfile = os.path.basename(args.eimage_file).replace('lsst_e', 'lsst_a')
-image_source.write_fits_file(outfile)
+persist = config['persistence']
+outfile = os.path.basename(args.eimage_file).replace(persist['eimage_prefix'],
+                                                     persist['raw_file_prefix'])
+
+image_source.write_fits_file(outfile, compress=persist['raw_file_compress'])

data/default_imsim_configs

@@ -12,6 +12,9 @@ bf_strength = 1.0
 [persistence]
 eimage_prefix = lsst_e_
 eimage_compress = True
+eimage_skip = True
+raw_file_prefix = lsst_a_
+raw_file_compress = True
 centroid_prefix = centroid_
 
 [cosmic_rays]

python/desc/imsim/ImageSimulator.py

@@ -20,6 +20,7 @@
 from .bleed_trails import apply_channel_bleeding
 from .skyModel import make_sky_model
 from .process_monitor import process_monitor
+from .camera_readout import ImageSource
 
 __all__ = ['ImageSimulator', 'compress_files']
 
@@ -97,8 +98,6 @@ def __init__(self, instcat, psf, numRows=None, config=None, seed=267,
         self._make_gs_interpreters(seed, sensor_list, file_id)
         self.log_level = log_level
         self.logger = get_logger(self.log_level, name='ImageSimulator')
-        if not self.gs_obj_arr:
-            self.logger.warning("No object entries in %s", instcat)
 
     def _make_gs_interpreters(self, seed, sensor_list, file_id):
         """
@@ -365,11 +364,14 @@ def __call__(self, gs_objects):
         if not os.path.isdir(outdir):
             os.makedirs(outdir)
         prefix = IMAGE_SIMULATOR.config['persistence']['eimage_prefix']
-        visit = str(IMAGE_SIMULATOR.obs_md.OpsimMetaData['obshistID'])
-        name_root = os.path.join(outdir, prefix) + visit
-        outfiles = gs_interpreter.writeImages(nameRoot=name_root)
-        if IMAGE_SIMULATOR.config['persistence']['eimage_compress']:
-            compress_files(outfiles)
+        obsHistID = str(IMAGE_SIMULATOR.obs_md.OpsimMetaData['obshistID'])
+        nameRoot = os.path.join(outdir, prefix) + obsHistID
+        if not IMAGE_SIMULATOR.config['persistence']['eimage_skip']:
+            outfiles = gs_interpreter.writeImages(nameRoot=nameRoot)
+            if IMAGE_SIMULATOR.config['persistence']['eimage_compress']:
+                compress_files(outfiles)
+        else:
+            self.write_raw_files(gs_interpreter.detectorImages)
 
         # Write out the centroid files if they were made.
         gs_interpreter.write_centroid_files()
@@ -400,6 +402,25 @@ def update_checkpoint_summary(self, gs_interpreter, num_objects):
                               gs_interpreter.nobj_checkpoint))
 
 
+    def write_raw_files(self, detector_images):
+        """
+        Write the raw files directly from galsim images.
+
+        Parameters
+        ----------
+        detector_images: dict
+            The dictionary attribute of the GalSimInterpreter object
+            used to fill the galsim images with eimage data.
+        """
+        persist = IMAGE_SIMULATOR.config['persistence']
+        prefix = persist['raw_file_prefix']
+        obsHistID = str(IMAGE_SIMULATOR.obs_md.OpsimMetaData['obshistID'])
+        nameRoot = os.path.join(IMAGE_SIMULATOR.outdir, prefix) + obsHistID
+        for name, gs_image in detector_images.items():
+            raw = ImageSource.create_from_galsim_image(gs_image)
+            outfile = '_'.join((nameRoot, name))
+            raw.write_fits_file(outfile, compress=persist['raw_file_compress'])
+
 def compress_files(file_list, remove_originals=True):
     """
     Use gzip to compress a list of files.

python/desc/imsim/camera_readout.py

@@ -21,6 +21,7 @@
 import scipy
 import astropy.io.fits as fits
 import astropy.time
+import galsim
 import lsst.afw.geom as afwGeom
 import lsst.afw.image as afwImage
 import lsst.utils as lsstUtils
@@ -30,6 +31,7 @@
     from lsst.sims.utils import getRotSkyPos
 from .camera_info import CameraInfo
 from .imSim import get_logger, get_config
+from .cosmic_rays import CosmicRays
 
 __all__ = ['ImageSource', 'set_itl_bboxes', 'set_e2v_bboxes',
            'set_phosim_bboxes', 'set_noao_keywords', 'cte_matrix']
@@ -57,18 +59,21 @@ class ImageSource(object):
         Object containing the readout properties of the sensors in the
         focal plane, provided by obs_lsstCam.ImsimMapper().camera.
     '''
-    def __init__(self, image_array, exptime, sensor_id, logger=None):
+    def __init__(self, image_array, exptime, sensor_id, visit=42, logger=None):
         """
         Class constructor.
 
         Parameters
         ----------
-        image_array : np.array
+        image_array: np.array
             A numpy array containing the pixel data for an eimage.
-        exptime : float
+        exptime: float
             The exposure time of the image in seconds.
-        sensor_id : str
+        sensor_id: str
             The raft and sensor identifier, e.g., 'R22_S11'.
+        visit: int [42]
+            Visit number to be used with the sensor_id for generating
+            the random seed for the dark current and read noise.
         logger: logging.Logger [None]
             logging.Logger object to use. If None, then a logger with level
             INFO will be used.
@@ -79,6 +84,7 @@ def __init__(self, image_array, exptime, sensor_id, logger=None):
 
         self.exptime = exptime
         self.sensor_id = sensor_id
+        self.visit = visit
 
         self.camera_info = CameraInfo()
 
@@ -96,6 +102,30 @@ def __init__(self, image_array, exptime, sensor_id, logger=None):
     def __del__(self):
         self.eimage.close()
 
+    @staticmethod
+    def create_from_galsim_image(gs_image, logger=None):
+        """
+        Create an ImageSource object from a galsim Image object
+        produced by a GalSimInterpreter.
+        """
+        # Extract the exptime and sensor_id from the
+        # sims_GalSimInterface.GalSim_afw_TanSimWCS object.
+        wcs = gs_image.wcs
+        exptime = wcs.photParams.exptime*wcs.photParams.nexp
+        sensor_id = 'R{}{}_S{}{}'\
+                    .format(*[x for x in wcs.detectorName if x.isdigit()])
+
+        raw_image = ImageSource(gs_image.array, exptime, sensor_id,
+                                visit=wcs.fitsHeader.getScalar('OBSID'))
+        raw_image.eimage_data = gs_image.array
+
+        # Add the keywords from wcs.fitsHeader to the raw_image.eimage
+        # attribute so that they are propagated to the raw file.
+        for name in wcs.fitsHeader.names():
+            raw_image.eimage[0].header[name] = wcs.fitsHeader.getScalar(name)
+        raw_image._read_pointing_info(None)
+        return raw_image
+
     @staticmethod
     def create_from_eimage(eimage_file, sensor_id=None, opsim_db=None,
                            logger=None):
@@ -132,6 +162,7 @@ def create_from_eimage(eimage_file, sensor_id=None, opsim_db=None,
             sensor_id = eimage[0].header['CHIPID']
 
         image_source = ImageSource(eimage[0].data, exptime, sensor_id,
+                                   visit=eimage[0].header['OBSID'],
                                    logger=logger)
         image_source.eimage = eimage
         image_source.eimage_data = eimage[0].data
@@ -149,19 +180,18 @@ def _read_pointing_info(self, opsim_db):
                 raise RuntimeError("eimage file does not have pointing info. "
                                    "Need an opsim db file.")
         # Read from the opsim db.
-        visit = self.eimage[0].header['OBSID']
         # We need an ObservationMetaData object to use the getRotSkyPos
         # function.
         obs_gen = ObservationMetaDataGenerator(database=opsim_db,
                                                driver="sqlite")
-        obs_md = obs_gen.getObservationMetaData(obsHistID=visit,
+        obs_md = obs_gen.getObservationMetaData(obsHistID=self.visit,
                                                 boundType='circle',
                                                 boundLength=0)[0]
         # Extract pointing info from opsim db for desired visit.
         conn = sqlite3.connect(opsim_db)
         query = """select descDitheredRA, descDitheredDec,
         descDitheredRotTelPos from summary where
-        obshistid={}""".format(visit)
+        obshistid={}""".format(self.visit)
         curs = conn.execute(query)
         ra, dec, rottelpos = [np.degrees(x) for x in curs][0]
         conn.close()
@@ -170,7 +200,6 @@ def _read_pointing_info(self, opsim_db):
         obs_md.pointingDec = dec
         self.rotangle = getRotSkyPos(ra, dec, obs_md, rottelpos)
 
-
     def get_amp_image(self, amp_info_record, imageFactory=afwImage.ImageI):
         """
         Return an amplifier afwImage.Image object with electronics
@@ -254,18 +283,23 @@ def _make_amp_image(self, amp_name):
             data = data[:, ::-1]
         if amp_info.getRawFlipY():
             data = data[::-1, :]
-
         imaging_segment.getArray()[:] = data
-        full_arr = full_segment.getArray()
 
         # Add dark current.
         dark_current = config['electronics_readout']['dark_current']
-        full_arr += np.random.poisson(dark_current*self.exptime,
-                                      size=full_arr.shape)
+        imaging_arr = imaging_segment.getArray()
+        # Generate a seed from the visit number and sensor_id so that
+        # the dark current noise is deterministic.
+        seed = CosmicRays.generate_seed(self.visit, self.sensor_id)
+        rng = galsim.PoissonDeviate(seed, dark_current*self.exptime)
+        dc_data = np.zeros(np.prod(imaging_arr.shape))
+        rng.generate(dc_data)
+        imaging_arr += dc_data.reshape(imaging_arr.shape)
 
         # Add defects.
 
         # Apply CTE.
+        full_arr = full_segment.getArray()
         pcti = config['electronics_readout']['pcti']
         pcte_matrix = cte_matrix(full_arr.shape[0], pcti)
         for col in range(0, full_arr.shape[1]):
@@ -293,8 +327,14 @@ def _add_read_noise_and_bias(self, amp_name):
         """
         amp_info = self.camera_info.get_amp_info(amp_name)
         full_arr = self.amp_images[amp_name].getArray()
-        full_arr += np.random.normal(scale=amp_info.getReadNoise(),
-                                     size=full_arr.shape)
+
+        # Generate a seed from the visit number and sensor_id so that
+        # the read noise is deterministic.
+        seed = CosmicRays.generate_seed(self.visit, self.sensor_id)
+        rng = galsim.GaussianDeviate(seed, amp_info.getReadNoise())
+        rn_data = np.zeros(np.prod(full_arr.shape))
+        rng.generate(rn_data)
+        full_arr += rn_data.reshape(full_arr.shape)
         full_arr += config['electronics_readout']['bias_level']
 
     def _apply_crosstalk(self):
@@ -313,28 +353,36 @@ def _apply_crosstalk(self):
             self.amp_images[amp_name].getArray()[:, :] \
                 = sum([x*y for x, y in zip(imarrs, xtalk_row)])
 
-    def get_amplifier_hdu(self, amp_name):
+    def get_amplifier_hdu(self, amp_name, compress=True):
         """
         Get an astropy.io.fits.HDU for the specified amplifier.
 
         Parameters
         ----------
-        amp_name : str
+        amp_name: str
             The amplifier name, e.g., "R22_S11_C00".
+        compress: bool [True]
+            Use RICE_1 compression.
 
         Returns
         -------
         astropy.io.fits.ImageHDU
             Image HDU with the pixel data and header keywords
             appropriate for the requested sensor segment.
         """
-        hdu = fits.ImageHDU(data=self.amp_images[amp_name].getArray().astype(np.int32))
+        data = self.amp_images[amp_name].getArray().astype(np.int32)
+        if compress:
+            hdu = fits.CompImageHDU(data=data, compression_type='RICE_1')
+        else:
+            hdu = fits.ImageHDU(data=data)
         hdr = hdu.header
         amp_info = self.camera_info.get_amp_info(amp_name)
         # Copy keywords from eimage primary header.
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
             for key in self.eimage[0].header.keys():
+                if key in ('BITPIX', 'NAXIS'):
+                    continue
                 try:
                     hdr[key] = self.eimage[0].header[key]
                 except ValueError:
@@ -390,21 +438,25 @@ def write_amplifier_image(self, amp_name, outfile, overwrite=True):
         output.writeto(outfile, overwrite=overwrite)
 
     def write_fits_file(self, outfile, overwrite=True, run_number=None,
-                        lsst_num='LCA-11021_RTM-000'):
+                        lsst_num='LCA-11021_RTM-000', compress=True):
         """
         Write the processed eimage data as a multi-extension FITS file.
 
         Parameters
         ----------
-        outfile : str
+        outfile: str
             Name of the output FITS file.
-        overwrite : bool, optional
-            Flag whether to overwrite an existing output file.  Default: True.
+        overwrite: bool [True]
+            Flag whether to overwrite an existing output file.
+        run_number: int [None]
+            Run number.  If None, then the visit number is used.
+        compress: bool [True]
+            Use RICE_1 compression for each image HDU.
         """
         output = fits.HDUList(fits.PrimaryHDU())
         output[0].header = self.eimage[0].header
         if run_number is None:
-            run_number = output[0].header['OBSID']
+            run_number = self.visit
         output[0].header['RUNNUM'] = str(run_number)
         output[0].header['DARKTIME'] = output[0].header['EXPTIME']
         mjd_obs = astropy.time.Time(output[0].header['MJD-OBS'], format='mjd')
@@ -428,7 +480,7 @@ def write_fits_file(self, outfile, overwrite=True, run_number=None,
         seg_ids = '10 11 12 13 14 15 16 17 07 06 05 04 03 02 01 00'.split()
         for seg_id in seg_ids:
             amp_name = '_C'.join((self.sensor_id, seg_id))
-            output.append(self.get_amplifier_hdu(amp_name))
+            output.append(self.get_amplifier_hdu(amp_name, compress=compress))
             output[-1].header['EXTNAME'] = 'Segment%s' % seg_id
         output.writeto(outfile, overwrite=overwrite)
 

tests/test_camera_readout.py

@@ -20,7 +20,12 @@ class ImageSourceTestCase(unittest.TestCase):
     image_source \
         = desc.imsim.ImageSource.create_from_eimage(eimage_file, 'R22_S11')
     def setUp(self):
-        pass
+        imsim_dir = lsstUtils.getPackageDir('imsim')
+        self.eimage_file = os.path.join(imsim_dir, 'tests', 'data',
+                                        'lsst_e_161899_R22_S11_r.fits.gz')
+        self.image_source \
+            = desc.imsim.ImageSource.create_from_eimage(self.eimage_file,
+                                                        'R22_S11')
 
     def tearDown(self):
         pass