Merge pull request #12 from spacetelescope/flatfield

Add support for flat field and gain reference files.

docs/index.rst

@@ -15,6 +15,7 @@ Contents
    romanisim/image
    romanisim/l1
    romanisim/l2
+   romanisim/refs
 
    romanisim/catalog
    romanisim/apt

docs/romanisim/refs.rst

@@ -0,0 +1,32 @@
+Reference files
+===============
+
+romanisim uses reference files from `CRDS <https://hst-crds.stsci.edu/static/users_guide/index.html>`_ in order to simulate realistic images.  The following kinds of reference files are used:
+
+* read noise
+* dark current
+* flat field
+* gain
+* distortion map
+
+The usage of these is mostly straightforward, but we provide here a few notes.
+
+Read Noise
+----------
+The random noise on reading a sample contributing to a ramp in an L1 image is scaled by the read noise reference file.
+
+Dark Current
+------------
+CRDS provides dark current images for each possible MA table, including the averaging of the dark current into resultants.  This simplifies subtraction from L1 images and allows effects beyond a simple Poisson sampling of dark current electrons in each read.  But it's unwieldy for a simulator because any effects beyond simple Poisson sampling of dark current electrons are not presently defined well enough to allow simulation.  So the simulator simply takes the last resultant in the dark current resultant image and scales it by the effective exposure time of that resultant to get a dark current rate.  This rate then goes into the idealized "counts" image which is then apportioned into the reads making up the resultants of an L1 image.
+
+Flat field
+----------
+Implementation of the flat field requires a little care due to the desire to support galsim's "photon shooting" rendering mode.  This mode does not create noise-free images but instead only simulates the number of photons that would be actually detected in a device.  We want to start by simulating the number of photons each pixel would record for a flat field of 1, and then sample that down by a fraction corresponding to the actual QE of each pixel.  That works fine supposing that the flat field is less than 1, but does not work for values of the flat field greater than 1.  So we instead do the initial galsim simulations for a larger exposure time than necessary, scaled by the maximum value of the flat field, and then sample down by :math:`\mathrm{flat}/\mathrm{maxflat}`.  That's all well and good as long as there aren't any spurious very large values in the flat field.  I haven't actually seen any such values yet and so haven't tried to address that case (e.g., by clipping them).
+
+Gain
+----
+Photons from the idealized "counts" image are scaled down to ADU before quantization during L1 creation, and then converted back to electrons before ramp fitting when making L2 images.
+
+Distortion map
+--------------
+World coordinate systems for WFI images are created by placing the telescope boresight at :math:`\mathrm{V2} = \mathrm{V3} = 0`, and then applying the distortion maps from CRDS to convert from V2V3 to pixels.

romanisim/catalog.py

@@ -64,7 +64,7 @@ def make_dummy_catalog(coord, radius=0.1, rng=None, seed=42, nobj=1000,
     objlist = []
     locs = util.random_points_in_cap(coord, radius, nobj, rng=rng)
     for i in range(nobj):
-        sky_pos = locs[i]
+        sky_pos = util.celestialcoord(locs[i])
         p = rng()
         # prescription follows galsim demo13.
         if p < 0.8:  # 80% of targets; faint galaxies
@@ -157,8 +157,8 @@ def make_dummy_table_catalog(coord, radius=0.1, rng=None, nobj=1000,
     hlr[star] = 0
 
     out = table.Table()
-    out['ra'] = [x.ra.deg for x in locs]
-    out['dec'] = [x.dec.deg for x in locs]
+    out['ra'] = locs.ra.to(u.deg).value
+    out['dec'] = locs.dec.to(u.deg).value
     out['type'] = types
     out['n'] = sersic_index
     out['half_light_radius'] = hlr

romanisim/image.py

@@ -51,8 +51,13 @@
 # I've at least caught these for the moment, but we should explore lower
 # folding_ratio for speed, as well as directly rendering into the image.
 
+# should we let people specify reference files?
+# what reference files do we use?
+# distortion, read noise, dark, gain, flat
+# these would each be optional arguments that would override
 
-def make_l2(resultants, ma_table, read_noise=None):
+
+def make_l2(resultants, ma_table, read_noise=None, gain=None, flat=None):
     """
     Simulate an image in a filter given resultants.
 
@@ -66,6 +71,8 @@ def make_l2(resultants, ma_table, read_noise=None):
         list of list of first read numbers and number of reads in each resultant
     read_noise : ndarray-like[n_pix, n_pix]
         read_noise image to use.  If None, use galsim.roman.read_noise.
+    flat : np.ndarray[float]
+        flat field to use
 
     Returns
     -------
@@ -80,18 +87,27 @@ def make_l2(resultants, ma_table, read_noise=None):
     if read_noise is None:
         read_noise = galsim.roman.read_noise
 
+    if gain is None:
+        gain = galsim.roman.gain
+
     from . import ramp
     rampfitter = ramp.RampFitInterpolator(ma_table)
     log.warning('Gain should be handled as something more interesting '
                 'than a single constant.')
-    ramppar, rampvar = rampfitter.fit_ramps(resultants/galsim.roman.gain,
-                                            read_noise)
+    ramppar, rampvar = rampfitter.fit_ramps(resultants*gain, read_noise)
     # could iterate if we wanted to improve the flux estimates
 
-    return (ramppar[..., 1],  # the slopes, ignoring the pedestals
-            rampvar[..., 0, 1, 1], # the read noise induced slope variance
-            rampvar[..., 1, 1, 1], # the poisson noise induced slope variance
-           )
+    slopes = ramppar[..., 1]
+    readvar = rampvar[..., 0, 1, 1]
+    poissonvar = rampvar[..., 1, 1, 1]
+
+    if flat is not None:
+        flat = np.clip(flat, 1e-9, np.inf)
+        slopes /= flat
+        readvar /= flat**2
+        poissonvar /= flat**2
+
+    return slopes, readvar, poissonvar
 
 
 def in_bounds(objlist, wcs, imbd, margin):
@@ -130,7 +146,7 @@ def simulate_counts(sca, targ_pos, date, objlist, filter_name,
                     exptime=None, rng=None, seed=None,
                     ignore_distant_sources=10, usecrds=True,
                     return_info=False, webbpsf=True,
-                    darkrate=None):
+                    darkrate=None, flat=None):
     """Simulate total counts in a single SCA.
 
     This gives the total counts in an idealized instrument with no systematics;
@@ -158,8 +174,10 @@ def simulate_counts(sca, targ_pos, date, objlist, filter_name,
         do not render sources more than this many pixels off edge of detector
     usecrds : bool
         use CRDS distortion map
-    darkrate : float or np.ndarray of float
+    darkrate : float or np.ndarray[float]
         dark rate image to use (electrons / s)
+    flat : float or np.ndarray[float]
+        flat field to use
 
     Returns
     -------
@@ -202,8 +220,22 @@ def simulate_counts(sca, targ_pos, date, objlist, filter_name,
     sky_image += roman.thermal_backgrounds[galsim_filter_name] * exptime
     imbd = full_image.bounds
     abflux = romanisim.bandpass.get_abflux(filter_name)
-    log.info('Adding sources to image...')
 
+    if flat is None:
+        flat = 1
+    # for some reason, galsim doesn't like multiplying an SED by 1, but it's
+    # okay with multiplying an SED by 1.0.
+    maxflat = float(np.max(flat))
+    if maxflat > 1.1:
+        log.warning('max(flat) > 1.1; this seems weird?!')
+    if maxflat > 2:
+        log.error('max(flat) > 2; this seems really weird?!')
+    # how to deal with the flat field?  We artificially inflate the
+    # exposure time of each source by maxflat when rendering.  And then we
+    # do a binomial sampling of the total number of photons obtained per pixel
+    # to figure out how many "should" have entered the pixel.
+
+    log.info('Adding sources to image...')
     nrender = 0
     final = None
     info = []
@@ -216,15 +248,15 @@ def simulate_counts(sca, targ_pos, date, objlist, filter_name,
             info.append(0)
             continue
         image_pos = galsim.PositionD(xpos[i], ypos[i])
-        final = galsim.Convolve(obj.profile * exptime, psf)
+        final = galsim.Convolve(obj.profile * exptime * maxflat, psf)
         if chromatic:
             stamp = final.drawImage(
                 bandpass, center=image_pos, wcs=imwcs.local(image_pos),
                 method='phot', rng=rng)
         else:
             if obj.flux is not None:
                 final = final.withFlux(
-                    obj.flux[filter_name] * abflux * exptime)
+                    obj.flux[filter_name] * abflux * exptime * maxflat)
                 try:
                     stamp = final.drawImage(center=image_pos,
                                             wcs=imwcs.local(image_pos))
@@ -249,12 +281,17 @@ def simulate_counts(sca, targ_pos, date, objlist, filter_name,
     if final is not None and not final.spectral:
         full_image.addNoise(poisson_noise)
 
+    if not np.all(flat == 1):
+        full_image.array[:, :] = np.random.binomial(
+            np.round(full_image.array).astype('i4'), flat/maxflat)
+
     full_image += sky_image
 
     if darkrate is None:
         darkrate = roman.dark_current
     dark_image += darkrate * exptime
     dark_image.addNoise(poisson_noise)
+    full_image += dark_image
 
     full_image.quantize()
 
@@ -322,18 +359,39 @@ def simulate(metadata, objlist,
     last_read = ma_table[-1][0] + ma_table[-1][1] - 1
     exptime = last_read * parameters.read_time
 
+    # TODO: replace this stanza with a function that looks at the metadata
+    # and keywords and returns a dictionary with all of the relevant reference
+    # data in numpy arrays.
+    # should query CRDS for any reference files not specified on the command
+    # line.
     if usecrds:
-        reffiles = crds.getreferences(all_metadata, reftypes=['readnoise', 'dark'],
-                                      observatory='roman')
+        reffiles = crds.getreferences(
+            all_metadata, observatory='roman',
+            reftypes=['readnoise', 'dark', 'gain'])
         read_noise = asdf.open(reffiles['readnoise'])['roman']['data']
         darkrate = asdf.open(reffiles['dark'])['roman']['data']
-        darkrate = darkrate[-1] / (np.mean(ma_table[-1])*parameters.read_time)
+        gain = asdf.open(reffiles['gain'])['roman']['data']
+        try:
+            reffiles = crds.getreferences(
+                all_metadata, observatory='roman',
+                reftypes=['flat'])
+            flat = asdf.open(reffiles['flat'])['roman']['data']
+        except crds.core.exceptions.CrdsLookupError:
+            log.warning('Could not find flat; using 1')
+            flat = 1
+        # convert the last dark resultant into a dark rate by dividing by the
+        # mean time in that resultant.
+        darkrate = darkrate[-1] / (
+            (ma_table[-1][0] + (ma_table[-1][1]-1)/2) * parameters.read_time)
         nborder = parameters.nborder
         read_noise = read_noise[nborder:-nborder, nborder:-nborder]
         darkrate = darkrate[nborder:-nborder, nborder:-nborder]
+        gain = gain[nborder:-nborder, nborder:-nborder]
     else:
         read_noise = galsim.roman.read_noise
         darkrate = galsim.roman.dark_current
+        gain = galsim.roman.gain
+        flat = 1
 
     out = dict()
 
@@ -349,15 +407,17 @@ def simulate(metadata, objlist,
     counts, simcatobj = simulate_counts(
         sca, coord, date, objlist, filter_name, rng=rng,
         exptime=exptime, usecrds=usecrds, darkrate=darkrate,
-        webbpsf=webbpsf)
+        webbpsf=webbpsf, flat=flat)
     if level == 0:
         return counts
     l1 = romanisim.l1.make_l1(
-        counts, ma_table_number, read_noise=read_noise, rng=rng, **kwargs)
+        counts, ma_table_number, read_noise=read_noise, rng=rng, gain=gain,
+        **kwargs)
     if level == 1:
         im = romanisim.l1.make_asdf(l1, metadata=all_metadata)
     elif level == 2:
-        slopeinfo = make_l2(l1, ma_table, read_noise=read_noise)
+        slopeinfo = make_l2(l1, ma_table, read_noise=read_noise,
+                            gain=gain, flat=flat)
         im = make_asdf(*slopeinfo, metadata=all_metadata)
     return im, simcatobj
 

romanisim/l1.py

@@ -164,7 +164,7 @@ def tij_to_pij(tij):
             tremaining -= (t - tlast)
             tlast = t
         pij.append(pi)
-    return pij
+    return np.clip(pij, 0, 1)
 
 
 def apportion_counts_to_resultants(counts, tij):
@@ -339,7 +339,8 @@ def ma_table_to_tij(ma_table_number):
 
 
 def make_l1(counts, ma_table_number,
-            read_noise=None, filepath=None, rng=None, seed=None):
+            read_noise=None, filepath=None, rng=None, seed=None,
+            gain=None):
     """Make an L1 image from a counts image.
 
     This apportions the total counts among the different resultants and adds
@@ -362,6 +363,8 @@ def make_l1(counts, ma_table_number,
         Random number generator to use
     seed : int
         Seed for populating RNG.  Only used if rng is None.
+    gain : float or np.ndarray[float]
+        Gain (electrons / count) for converting counts to electrons
 
     Returns
     -------
@@ -407,6 +410,6 @@ def make_l1(counts, ma_table_number,
     log.warning('We need to make sure we have the right units on the read '
                 'noise.')
 
-    resultants /= roman.gain
+    resultants /= gain
     resultants = np.round(resultants)
     return resultants

romanisim/ramp.py

@@ -347,7 +347,12 @@ def variances(self, flux, read_noise):
         fluxonreadvar = np.clip(
             fluxonreadvar, self.flux_on_readvar_pts[0],
             self.flux_on_readvar_pts[-1])
-        return self.var_interpolator(fluxonreadvar).astype('f4')*read_noise**2
+        var = self.var_interpolator(fluxonreadvar).astype('f4')
+        read_noise = np.array(read_noise)
+        read_noise = read_noise.reshape(
+            read_noise.shape + (1,)*(len(var.shape)-len(read_noise.shape)))
+        var *= read_noise**2
+        return var
 
     def fit_ramps(self, resultants, read_noise, fluxest=None):
         """Fit ramps for a set of resultants and their read noise.

romanisim/util.py

@@ -119,8 +119,7 @@ def random_points_in_cap(coord, radius, nobj, rng=None):
     dist = np.arccos(1-(1-np.cos(np.radians(radius)))*dist)
     c1 = SkyCoord(coord.ra.rad*u.rad, coord.dec.rad*u.rad, frame='icrs')
     c1 = c1.directional_offset_by(ang*u.rad, dist*u.rad)
-    sky_pos = [celestialcoord(x) for x in c1]
-    return sky_pos
+    return c1
 
 
 def flatten_dictionary(d):