Merge branch 'main' of https://github.com/spacetelescope/romanisim in…

…to SourceInjectionTestHotfix

README.md

@@ -3,7 +3,8 @@
 # romanisim: an image simulator for Roman
 
 romanisim is a Galsim-based simulator of imaging data from the Wide
-Field Instrument (WFI) on the Nancy Grace Roman Space Telescope.  It uses
+Field Instrument (WFI) on the Nancy Grace Roman Space Telescope
+(pronounced roman-eye-sim, stylized Roman I-Sim).  It uses
 [Galsim](https://galsim-developers.github.io/GalSim/_build/html/overview.html)
 to render astronomical scenes,
 [WebbPSF](https://galsim-developers.github.io/GalSim/_build/html/overview.html)

docs/index.rst

@@ -3,6 +3,9 @@ romanisim
 
 A Roman WFI image simulator based on galsim.
 
+We stylize the simulator Roman I-Sim and pronounce it roman - eye -
+sim; the package is romanisim.
+
 Contents
 ========
 

docs/romanisim/install.rst

@@ -60,4 +60,5 @@ before running the above, or choose a different directory for
 WebbPSF's data files.
 
 Some users report issues with the FFTW dependency of galsim on Mac Arm
-systems.  See galsim's installation page for hints there.
+systems.  See galsim's installation page for hints there.  In
+particular it may be helpful to install FFTW before galsim and romanisim.

romanisim/bandpass.py

@@ -5,12 +5,12 @@
 the nominal flat AB spectrum of 3631 Jy.  The ultimate source of this
 information is https://roman.gsfc.nasa.gov/science/WFI_technical.html .
 """
-import os
-import pkg_resources
+from importlib import resources
 from scipy import integrate
 from astropy import constants
 from astropy.table import Table
 from astropy import units as u
+import functools
 
 # 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
@@ -51,11 +51,17 @@ def read_gsfc_effarea(filename=None):
     astropy.table.Table
         table with effective areas for different Roman bandpasses.
     """
+
+    reader = functools.partial(Table.read, format='csv', delimiter=',',
+                               header_start=1, data_start=2)
     if filename is None:
-        dirname = pkg_resources.resource_filename('romanisim', 'data')
-        filename = os.path.join(dirname, 'Roman_effarea_20201130.csv')
-    return Table.read(filename, format='csv', delimiter=',', header_start=1,
-                      data_start=2)
+        ref = (resources.files('romanisim') / 'data'
+               / 'Roman_effarea_20201130.csv')
+        with resources.as_file(ref) as filename:
+            out = reader(filename)
+    else:
+        out = reader(filename)
+    return out
 
 
 def compute_abflux(effarea=None):

romanisim/image.py

@@ -21,6 +21,7 @@
 from . import parameters
 from . import util
 from . import nonlinearity
+from . import ramp
 import romanisim.l1
 import romanisim.bandpass
 import romanisim.psf
@@ -63,7 +64,7 @@
 # these would each be optional arguments that would override
 
 
-def make_l2(resultants, ma_table, read_noise=None, gain=None, flat=None,
+def make_l2(resultants, read_pattern, read_noise=None, gain=None, flat=None,
             linearity=None, dark=None, dq=None):
     """
     Simulate an image in a filter given resultants.
@@ -74,8 +75,8 @@ def make_l2(resultants, ma_table, read_noise=None, gain=None, flat=None,
     ----------
     resultants : np.ndarray[nresultants, nx, ny]
         resultants array
-    ma_table : list[list] (int)
-        list of list of first read numbers and number of reads in each resultant
+    read_pattern : list[list] (int)
+        list of list of indices of reads entering each resultant
     read_noise : np.ndarray[nx, ny] (float)
         read_noise image to use.  If None, use galsim.roman.read_noise.
     flat : np.ndarray[nx, ny] (float)
@@ -111,12 +112,11 @@ def make_l2(resultants, ma_table, read_noise=None, gain=None, flat=None,
     if dark is not None:
         resultants = resultants - dark
 
-    from . import ramp
     log.info('Fitting ramps.')
 
     # commented out code below is inverse-covariance ramp fitting
     # which doesn't presently support DQ information
-    # rampfitter = ramp.RampFitInterpolator(ma_table)
+    # rampfitter = ramp.RampFitInterpolator(read_pattern)
     # ramppar, rampvar = rampfitter.fit_ramps(resultants * gain,
     #                                         read_noise * gain)
 
@@ -129,7 +129,7 @@ def make_l2(resultants, ma_table, read_noise=None, gain=None, flat=None,
 
     ramppar, rampvar = ramp.fit_ramps_casertano(
         resultants * gain, dq & parameters.dq_do_not_use,
-        read_noise * gain, ma_table)
+        read_noise * gain, read_pattern)
 
     log.warning('The ramp fitter is unaware of noise from dark current because '
                 'it runs on dark-subtracted images.  We could consider adding '
@@ -509,11 +509,11 @@ def simulate_counts(metadata, objlist,
         catalog of simulated objects in image
     """
 
-    ma_table = parameters.ma_table[
+    read_pattern = parameters.read_pattern[
         metadata['exposure']['ma_table_number']]
 
     sca = int(metadata['instrument']['detector'][3:])
-    exptime = parameters.read_time * (ma_table[-1][0] + ma_table[-1][1] - 1)
+    exptime = parameters.read_time * read_pattern[-1][-1]
     if rng is None and seed is None:
         seed = 43
         log.warning(
@@ -622,9 +622,8 @@ def simulate(metadata, objlist,
     ma_table_number = image_mod.meta.exposure.ma_table_number
     filter_name = image_mod.meta.instrument.optical_element
 
-    ma_table = parameters.ma_table[ma_table_number]
-    exptime_tau = ((ma_table[-1][0] + (ma_table[-1][1] / 2))
-                   * parameters.read_time)
+    read_pattern = parameters.read_pattern[ma_table_number]
+    exptime_tau = np.mean(read_pattern[-1]) * 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
@@ -734,7 +733,7 @@ def simulate(metadata, objlist,
         im, extras = romanisim.l1.make_asdf(
             l1, dq=l1dq, metadata=image_mod.meta, persistence=persistence)
     elif level == 2:
-        slopeinfo = make_l2(l1, ma_table, read_noise=read_noise,
+        slopeinfo = make_l2(l1, read_pattern, read_noise=read_noise,
                             gain=gain, flat=flat, linearity=linearity,
                             dark=dark, dq=l1dq)
         l2dq = np.bitwise_or.reduce(l1dq, axis=0)

romanisim/l1.py

@@ -425,36 +425,24 @@ def make_asdf(resultants, dq=None, filepath=None, metadata=None, persistence=Non
     return out, extras
 
 
-def ma_table_to_tij(ma_table_number, read_time=None):
-    """Get the times of each read going into resultants for a MA table.
-
-    Currently only ma_table_number = 1 is supported, corresponding to a simple
-    fiducial high latitude imaging MA table.
-
-    This presently uses a hard-coded, somewhat inflexible MA table description
-    in the parameters file.  But that seems like an okay option given that the
-    current 'official' file is slated for redesign when the format is relaxed.
+def read_pattern_to_tij(read_pattern):
+    """Get the times of each read going into resultants for a read_pattern.
 
     Parameters
     ----------
-    ma_table_number : int or list[list]
-        if int, id of multiaccum table to use
-        otherwise a list of (first_read, n_reads) tuples going into resultants.
-    read_time : number
-        Frame read time to use.  If None, use romanisim.parameters.read_time.
+    read_pattern : int or list[list]
+        If int, id of ma_table to use.
+        Otherwise a list of lists giving the indices of the reads entering each
+        resultant.
 
     Returns
     -------
     list[list[float]]
         list of list of readout times for each read entering a resultant
     """
-    if read_time is None:
-        read_time = parameters.read_time
-    if isinstance(ma_table_number, int):
-        tab = parameters.ma_table[ma_table_number]
-    else:
-        tab = ma_table_number
-    tij = [read_time * np.arange(f, f + n) for (f, n) in tab]
+    if isinstance(read_pattern, int):
+        read_pattern = parameters.read_pattern[read_pattern]
+    tij = [parameters.read_time * np.array(x) for x in read_pattern]
     return tij
 
 
@@ -483,7 +471,7 @@ def add_ipc(resultants, ipc_kernel=None):
     return out
 
 
-def make_l1(counts, ma_table_number,
+def make_l1(counts, read_pattern,
             read_noise=None, rng=None, seed=None,
             gain=None, inv_linearity=None, crparam=None,
             persistence=None, tstart=None, saturation=None):
@@ -496,9 +484,9 @@ def make_l1(counts, ma_table_number,
     ----------
     counts : galsim.Image
         total counts delivered to each pixel
-    ma_table_number : int
-        multi accum table number indicating how reads are apportioned among
-        resultants
+    read_pattern : int or list[list]
+        MA table number or list of lists giving indices of reads entering each
+        resultant.
     read_noise : np.ndarray[nx, ny] (float) or float
         Read noise entering into each read
     rng : galsim.BaseDeviate
@@ -526,7 +514,7 @@ def make_l1(counts, ma_table_number,
         DQ array marking saturated pixels and cosmic rays
     """
 
-    tij = ma_table_to_tij(ma_table_number)
+    tij = read_pattern_to_tij(read_pattern)
     log.info('Apportioning counts to resultants...')
     resultants, dq = apportion_counts_to_resultants(
         counts.array, tij, inv_linearity=inv_linearity, crparam=crparam,

romanisim/parameters.py

@@ -5,6 +5,22 @@
 from astropy.time import Time
 from astropy import units as u
 
+read_pattern = {1: [[1 + x for x in range(8)],
+                    [9 + x for x in range(8)],
+                    [17 + x for x in range(8)],
+                    [25 + x for x in range(8)],
+                    [33 + x for x in range(8)],
+                    [41 + x for x in range(8)]],
+                2: [[1 + x for x in range(5)],
+                    [6 + x for x in range(8)],
+                    [14],
+                    [15 + x for x in range(9)],
+                    [24 + x for x in range(25)]],
+                3: [[1 + x for x in range(25)],
+                    [26 + x for x in range(8)],
+                    [34],
+                    [35 + x for x in range(14)]]
+                }
 
 default_parameters_dictionary = {
     'instrument': {'name': 'WFI',
@@ -14,6 +30,7 @@
     'exposure': {'start_time': Time('2026-01-01T00:00:00'),
                  'type': 'WFI_IMAGE',
                  'ma_table_number': 1,
+                 'read_pattern': read_pattern[1],
                  },
     'pointing': {'ra_v1': 270.0,
                  'dec_v1': 66.0,
@@ -45,10 +62,6 @@
 
 nborder = 4  # number of border pixels used for reference pixels.
 
-ma_table = {1: [[1, 8], [9, 8], [17, 8], [25, 8], [33, 8], [41, 8]],
-            2: [[1, 5], [6, 8], [14, 1], [15, 9], [24, 25]],
-            3: [[1, 25], [26, 8], [34, 1], [35, 14]]}
-
 read_time = 3.04
 
 # from draft wfisim documentation