add script to create flats

bin.src/make_flats.py

@@ -0,0 +1,51 @@
+#!/usr/bin/env python
+"""
+Script to create imSim flats.  This includes electorstatic sensor effects
+such as treerings and brighter-fatter.
+"""
+import argparse
+import galsim
+import lsst.afw.cameraGeom as cameraGeom
+from lsst.sims.GalSimInterface import LSSTCameraWrapper, make_galsim_detector
+import desc.imsim
+
+parser = argparse.ArgumentParser("imSim flat production script.")
+parser.add_argument('instcat', type=str, help="instance catalog file")
+parser.add_argument('--sensors', type=str, default=None,
+                    help="sensors to simulate")
+parser.add_argument('--seed', type=int, default=42, help="random seed")
+parser.add_argument('--counts_total', type=int, default=8e4,
+                    help="total # electrons/pixel in flat")
+parser.add_argument('--counts_per_iter', type=int, default=4e3,
+                    help="# electrons/pixel per iteration")
+parser.add_argument('--log_level', type=str,
+                    choices=['DEBUG', 'INFO', 'WARN', 'ERROR', 'CRITICAL'],
+                    default='INFO', help='Logging level. Default: INFO')
+
+args = parser.parse_args()
+
+obs_md, phot_params, _ = desc.imsim.parsePhoSimInstanceFile(args.instcat,
+                                                            numRows=30)
+sensor_list = args.sensors.split('^') if args.sensors is not None \
+    else args.sensors
+rng = galsim.UniformDeviate(args.seed)
+niter = int(args.counts_total/args.counter_per_iter + 0.5)
+counts_per_iter = args.counts_total/niter
+logger = desc.imsim.get_logger(args.log_level, name='make_flats')
+config = desc.imsim.read_config()
+
+visit = obs_md.OpsimMetaData['obshistID']
+
+camera_wrapper = LSSTCameraWrapper()
+
+for det in camera_wrapper.camera:
+    det_name = det.getName()
+    logger.info("processing %s", det_name)
+    if (det.getType() != cameraGeom.SCIENCE or
+            (args.sensors is not None and det_name not in sensor_list)):
+        continue
+    gs_det = make_galsim_detector(camera_wrapper, det_name, phot_params, obs_md)
+    desc.imsim.add_treering_info([gs_det])
+    my_flat = desc.imsim.make_flat(gs_det, counts_per_iter, niter, rng)
+    ccd_id = "R{}_S{}".format(det_name[2:5:2], det_name[8:11:2])
+    my_flat.write('flat_{}_{}_{}.fits'.format(visit, ccd_id, obs_md.bandpass))

python/desc/imsim/__init__.py

@@ -18,3 +18,4 @@
 from .bleed_trails import *
 from .process_monitor import *
 from .instcat_tools import *
+from .flats import *

python/desc/imsim/flats.py

@@ -0,0 +1,77 @@
+"""
+imSim module to make flats.
+"""
+import galsim
+
+__all__ = ['make_flat']
+
+def make_flat(gs_det, counts_per_iter, niter, rng, buf=2):
+    """
+    Create a flat by successively adding lower level flat sky images.
+    This is based on
+    https://github.com/GalSim-developers/GalSim/blob/releases/2.0/devel/lsst/treering_flat.py
+
+    Full LSST CCDs are assembled one amp at a time to limit the memory
+    used by the galsim.SiliconSensor model.
+
+    Parameters
+    ----------
+    gs_det: GalSimDetector
+        The detector in the LSST focalplane to use.  This object
+        contains the CCD pixel geometry, WCS, and treering info.
+    counts_per_iter: int
+        Roughly the number of electrons/pixel to add at each iteration.
+    niter: int
+        Number of iterations. Final flat will have niter*counts_per_iter
+        electrons/pixel.
+    rng: galsim.BaseDeviate
+        Random number generator.
+    buf: int [2]
+        Pixel buffer around each to account for charge redistribution
+        across pixel boundaries.
+
+    Returns
+    -------
+    galsim.ImageF
+    """
+    ncol = gs_det.xMaxPix - gs_det.xMinPix + 1
+    nrow = gs_det.yMaxPix - gs_det.yMinPix + 1
+    flat = galsim.ImageF(ncol, nrow, wcs=gs_det.wcs)
+    sensor = galsim.SiliconSensor(rng=rng,
+                                  treering_center=gs_det.tree_rings.center,
+                                  treering_func=gs_det.tree_rings.func,
+                                  transpose=True)
+
+    # Create a noise-free base image to add at each iteration.
+    base_image = galsim.ImageF(ncol, nrow, wcs=gs_det.wcs)
+    base_image.wcs.makeSkyImage(base_image, sky_level=1.)
+    mean_pixel_area = base_image.array.mean()
+    sky_level_per_iter = counts_per_iter/mean_pixel_area
+    base_image *= sky_level_per_iter
+
+    noise = galsim.PoissonNoise(rng)
+
+    # Build up the full CCD by amplifier segment to limit the memory
+    # used by the silicon model.
+    nx, ny = 2, 8
+    dx = ncol//nx
+    dy = nrow//ny
+    for i in range(nx):
+        xmin = i*dx + 1
+        xmax = (i + 1)*dx
+        for j in range(ny):
+            ymin = j*dy + 1
+            ymax = (j + 1)*dy
+            temp_image = galsim.ImageF(ncol, nrow, wcs=gs_det.wcs)
+            bounds = (galsim.BoundsI(xmin-buf, xmax+buf, ymin-buf, ymax+buf)
+                      & temp_image.bounds)
+            temp_amp = temp_image[bounds]
+            for _ in range(niter):
+                temp = sensor.calculate_pixel_areas(temp_amp)
+                temp /= temp.mean()
+                temp *= base_image[bounds]
+                temp.addNoise(noise)
+                temp_amp += temp
+            amp_bounds = galsim.BoundsI(xmin, xmax, ymin, ymax)
+            flat[amp_bounds] += temp_image[amp_bounds]
+    return flat