Add romanisim-make-l3 image; refactor L3 code to support (#145)

* Initial L3 refactoring commit.
* Add romanisim-make-l3.
* Add L3 documentation.
* Add L3 metadata.

docs/index.rst

@@ -19,6 +19,7 @@ Contents
    romanisim/image
    romanisim/l1
    romanisim/l2
+   romanisim/l3
    romanisim/refs
 
    romanisim/catalog

docs/romanisim/image.rst

@@ -14,4 +14,19 @@ The parsed metadata is used to make a ``counts`` image that is an idealized imag
 These idealized count images are then used to either make a level 2 image or a level 1 image, which are intended to include
 the effects of these complications.  The construction of L1 images is described :doc:`here </romanisim/l1>`, and the construction of L2 images is described :doc:`here </romanisim/l2>`.
 
+L2 source injection
+-------------------
+
+We support injecting sources into L2 files.  L2 source injection follows many of the same steps as L1 creation, with some short cuts.  The simulation creates an idealized image of the total number of counts deposited into each pixel, using the provided catalog, exposure time, filter, WCS, and PSF from the provided image.  A virtual ramp is generated for the existing L2 file by evenly apportioning the L2 flux among the resultants.  Additional counts from the injected source are added to each resultant in the virtual ramp using the same code as for the L1 simulations.  The resulting virtual ramp is fit to get the new flux per pixel, and we replace the values in the original L2 model with the new slope measurements.
+
+This makes some shortcuts:
+
+* The L2 files don't include information about which pixels contain CR hits.  Sources injected into pixels with CR hits get re-fit as if there were no CR.
+* Non-linearity is not included; the ramp is refit without any non-linearity.
+* The ramp is refit without persistence, but any persistence which was included in the initial ramp slope is still included.
+
+The function :meth:`romanisim.image.inject_sources_into_l2` is the intended entry point for L2 source injection; see its documentation for more details about how to add sources to an L2 image.
+
+
+
 .. automodapi:: romanisim.image

docs/romanisim/l2.rst

@@ -16,4 +16,16 @@ This is a fairly faithful representation of how level two image construction wor
 * We have a simplistic saturation treatment, just clipping resultants that exceed
   the saturation level to the saturation level and throwing a flag.
 
+
+L2 source injection
+===================
+We support L2 source injection through :meth:`romanisim.image.inject_sources_into_l2`.  This routine is intended for cases when existing L2 images are available, and a small number of sources need to be added on top of those images.  The L2 source injection is intended to be a fairly faithful simulation of the process.  It includes the following steps:
+* The model WCS is used to generate the pixel locations for sources.
+* The number of additional electrons in each pixel is simulated.
+* A virtual ramp is created that evenly apportions the original L2 flux along the ramp.
+* The new counts are apportioned randomly along the ramp following a Poisson distribution in the same way as for the usual L1 simulations.
+* The new ramp is then fit using the ramp fitting code, and the new fit replaces the old fit in the L2 file.
+
+Note that this procedure does not include any new cosmic rays, or any special treatment of pixels that had cosmic rays in them and now have additional source flux.  It is hard to do this right since we do not know which pixel the cosmic ray landed in on the basis of the L2 file.  If you need more accuracy than this you probably need to directly inject sources into the L1 file.
+
 .. automodapi:: romanisim.ramp

docs/romanisim/l3.rst

@@ -0,0 +1,59 @@
+L3 images
+=========
+
+We also simulate L3 images, also known as coadded images or mosaics; these terms are used equivalently here.  In contrast to L1 and L2 images, L3 images do not do a full simulation of the individual reads entering each resultant and the following ramp fitting.  Accordingly, they also do not include most systematic effects---for example, cosmic rays and non-linearity are not simulated.  Instead, the L3 simulations conceptually are more closely related to the initial idealized :doc:`count images </romanisim/image>` than they are to L1 or L2 images.  They differ from these idealized images in that we attempt to change the PSF in order to account for non-native pixel scales that are possible in coadded images.
+
+The simulation of L3 images is simplified relative to true L3 images.  First, as above, we do not do the full L2 simulations.  Second, we do not have a notion of a series of images that are being coadded---we are instead just saying that one could produce a mosaic that has a particular pixel scale, effective read noise, and exposure time.  We do not attempt to simulate the drizzle process used for typical L3 images in Roman, and accordingly do not include any drizzle-induced covariances among the pixels.  Similarly, because we don't really have a notion of a series of input L2 images, the L3 PSF is not a combination of the PSF of the input images, but is instead taken as the PSF evaluated at the center of the WFI02 detector, and then additionally convolved depending on the output pixel grid scale.
+
+The L3 PSF aims to include the convolution with the native scale that Roman necessarily includes, but to also support oversampled images where the L3 pixel scale can be significantly finer than the native pixel scale.  In these cases we include an additional square convolution to "make up the difference" between the L3 pixel scale and the L2 pixel scale; see the :meth:`romanisim.l3.l3_psf` method for more details.
+
+The ``romanisim-make-l3`` command-line interface can be used to make L3 files.  Invocations look like::
+
+    romanisim-make-l3 --bandpass F158 --radec 270 66 --npix 2000 --pixscalefrac 0.9 --exptime 1000 l3.asdf gaia-270-66-2027-06-01.ecsv
+
+which would make an F158 image with 2000x2000 pixels, slightly oversampled, corresponding to 1000 seconds, using the given catalog and output filename.
+
+There are a number of possible input arguments::
+
+    romanisim-make-l3 -h
+    usage: romanisim-make-l3 [-h] [--bandpass BANDPASS] [--config CONFIG] [--date DATE] [--nobj NOBJ] [--radec RADEC RADEC] [--npix NPIX]
+			     [--pixscalefrac PIXSCALEFRAC] [--exptime EXPTIME] [--effreadnoise EFFREADNOISE] [--nexposures NEXPOSURES]
+			     [--rng_seed RNG_SEED]
+			     filename catalog
+
+    Make an L3 image.
+
+    positional arguments:
+      filename              output image (asdf)
+      catalog               input catalog (ecsv)
+
+    options:
+      -h, --help            show this help message and exit
+      --bandpass BANDPASS   bandpass to simulate (default: F087)
+      --config CONFIG       input parameter override file (yaml) (default: None)
+      --date DATE           UTC Date and Time of observation to simulate in ISOT format. (default: None)
+      --radec RADEC RADEC   ra and dec (deg) (default: None)
+      --npix NPIX           number of pixels across image (default: 4000)
+      --pixscalefrac PIXSCALEFRAC
+			    pixel scale as fraction of original (default: 1.0)
+      --exptime EXPTIME     total exposure time on field; roughly time per exposure times number of exposures (default: 100.0)
+      --effreadnoise EFFREADNOISE
+			    effective readnoise per pixel in MJy/sr. If None, a pessimistic estimate is computed. (default: None)
+      --nexposures NEXPOSURES
+			    number of exposures on field. Used only to compute the effective read noise. (default: 1)
+      --rng_seed RNG_SEED
+
+    EXAMPLE: romanisim-make-l3 output_image.asdf catalog.ecsv
+
+The arguments control things like the filter used, the date of observation (used to estimate the sky background), the right ascension and declination of the center of the mosaic, the number of pixels in the mosaic, and the degree of oversampling of the output mosaic.  Parameters also control the effective depth of the mosaic, through the exposure time, effective read noise, and number of exposures.
+
+L3 source injection
+-------------------
+
+We support injecting sources into L3 files.  L3 source injection is very similar to L3 image creation, except we do not need to invent effective read noises or exposure times or WCSes.  Instead we only change the Poisson noise, since we are only adding additional new sources.
+
+Conceptually we just add more sources to the image in the same way as in L3 image creation.  The primary challenge is picking the correct exposure time so that the Poisson noise in the simulation is correct; we (optionally) simulate individual photons and so the conversion between calibrated units and electrons is still relevant in the L3 files.  We do this by using the relationship between the calibrated flux and corresponding Poisson variance to figure out the effective calibrated-flux-to-electron relationship, and then simulate the appropriate number of electrons and convert back.
+
+The function :meth:`romanisim.l3.inject_sources_into_l3` is the intended entry point for L3 source injection; see its documentation for more details about how to add sources to an L3 image.
+
+.. automodapi:: romanisim.l3

romanisim/bandpass.py

@@ -12,6 +12,7 @@
 from astropy.table import Table
 from astropy import units as u
 import functools
+from romanisim import parameters
 
 # to go from calibrated fluxes in maggies to counts in the Roman bands
 # we need to multiply by a constant determined by the AB magnitude
@@ -181,3 +182,26 @@ def compute_count_rate(flux, bandpass, filename=None, effarea=None, wavedist=Non
     # per second
 
     return zpflux
+
+
+def etomjysr(bandpass):
+    """Compute factor converting e/s/pix to MJy/sr.
+
+    Assumes a pixel scale of 0.11" (romanisim.parameters.pixel_scale)
+
+    Parameters
+    ----------
+    bandpass : str
+        the name of the bandpass
+
+    Returns
+    -------
+    float
+        the factor F such that MJy / sr = F * DN/s
+    """
+
+    abflux = get_abflux(bandpass)  # e/s corresponding to 3631 Jy
+    srpix = ((parameters.pixel_scale * u.arcsec) ** 2).to(u.sr).value
+    mjysr = 1 / abflux * 3631 / 10 ** 6 / srpix  # should be ~0.3
+    return mjysr
+

romanisim/image.py

@@ -205,7 +205,7 @@ def trim_objlist(objlist, image):
 
 
 def add_objects_to_image(image, objlist, xpos, ypos, psf,
-                         flux_to_counts_factor, convtimes=None,
+                         flux_to_counts_factor, outputunit_to_electrons=None,
                          bandpass=None, filter_name=None, add_noise=False,
                          rng=None, seed=None):
     """Add sources to an image.
@@ -227,8 +227,9 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf,
         physical fluxes in objlist (whether in profile SEDs or flux arrays)
         should be multiplied by this factor to convert to total counts in the
         image
-    convtimes: array_like
-        Exposure times with unit scaling to convert to output rate units
+    outputunit_to_electrons : array_like
+        One output image unit corresponds to this many electrons.  If None,
+        leave as electrons.
     bandpass : galsim.Bandpass
         bandpass in which image is being rendered.  This is used only in cases
         where chromatic profiles & PSFs are being used.
@@ -287,17 +288,17 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf,
             stamp = final.drawImage(
                 bandpass, center=image_pos, wcs=pwcs,
                 method='phot', rng=rng)
-            if add_noise:
-                stamp.addNoise(galsim.PoissonNoise(rng))
         else:
             try:
                 stamp = final.drawImage(center=image_pos,
                                         wcs=pwcs)
+                if add_noise:
+                    stamp.addNoise(galsim.PoissonNoise(rng))
             except galsim.GalSimFFTSizeError:
                 log.warning(f'Skipping source {i} due to too '
                             f'large FFT needed for desired accuracy.')
-        if convtimes is not None:
-            stamp /= convtimes[i]
+        if outputunit_to_electrons is not None:
+            stamp /= outputunit_to_electrons[i]
 
         bounds = stamp.bounds & image.bounds
         if bounds.area() > 0:
@@ -599,7 +600,8 @@ def gather_reference_data(image_mod, usecrds=False):
     """
 
     reffiles = {k: v for k, v in parameters.reference_data.items()}
-    reffiles.pop('photom')
+    if 'photom' in reffiles:
+        reffiles.pop('photom')
 
     out = dict(**reffiles)
     if usecrds:
@@ -917,7 +919,7 @@ def inject_sources_into_l2(model, cat, x=None, y=None, psf=None, rng=None,
     This routine allows sources to be injected onto an existing L2 image.
     Source injection into an L2 image relies on knowing the objects'
     x and y locations, the PSF, and the image gain; if these are not provided,
-    reasonable defaults are generated.
+    reasonable defaults are generated from the input model.
 
     The simulation proceeds by (optionally) using the model WCS to generate the
     x & y locations, grabbing the gain from
@@ -1036,7 +1038,7 @@ def inject_sources_into_l2(model, cat, x=None, y=None, psf=None, rng=None,
         newramp[:, m], read_pattern,
         gain=gain, flat=1, darkrate=0)
 
-    res = model.copy()
+    res = copy.deepcopy(model)
     res.data[m] = newimage
     res.var_rnoise[m] = readvar
     res.var_poisson[m] = poissonvar