romanisim¶
A Roman WFI image simulator based on galsim.
Contents¶
Overview¶
romanisim simulates Roman Wide-Field Imager images of astronomical scenes described by catalogs of sources. The simulation includes:
convolution of the sources by the appropriate PSF for each detector
realistic world coordinate system
sky background
level 1 image support (3D image stack of up-the-ramp samples)
level 2 image support (2D image after calibration & ramp fitting)
point sources and analytic galaxy profiles
expected system throughput
dark current
read noise
inter-pixel capacitance
non-linearity
reciprocity failure
The simulator is based on galsim and most of these features directly invoke the equivalents in the galsim.roman package. The chief additions of this package on top of the galsim.roman implementation are using “official” PSF, WCS, and reference files from the Roman CRDS (not yet public!) and webbpsf. This package also implements WFI up-the-ramp sampled and averaged images like those that will be downlinked from the telescope, and the official Roman WFI file format (asdf).
The best way to interact with romanisim is to make an image. Running
romanisim-make-image out.asdf
will make a test image in the file out.asdf. Naturally, usually one has a particular astronomical scene in mind, and one can’t really simulate a scene without knowing where the telescope is pointing and when the observation is being made. A more complete invocation would be
romanisim-make-image –catalog input.ecsv –radec 270 66 –bandpass F087 –sca 7 –date 2026 1 1 –level 1 out.asdf
where input.ecsv includes a list of sources in the scene, the telescope boresight is pointing to (r, d) = (270, 66), the desired bandpass is F087, the sensor is WFI07, the date is Jan 1, 2026, and a level 1 image (3D cube of samples up the ramp) is requested.
Features not included so far:
any pedestal/frame 0 features
non-linear dark features
realistic treatment of the noise properties in L2 images, using information about how we go from L1 images & resultants to rate images.
Installation¶
First, install galsim, paying special attention to the dependence on FFTW. See the documentation here.
Then
pip install romanisim
and you should be set!
Making images¶
Ultimately, romanisim builds image by feeding source profiles, world coordinate system objects, and point spread functions to galsim. The image and l1 modules currently implement this functionality.
The image module is responsible for translating a metadata object that specifies everything about the conditions of the observation into objects that the simulation can understand. The metadata object follows the metadata that real WFI images will include; see here <https://roman-pipeline.readthedocs.io/en/latest/roman/datamodels/metadata.html#metadata> for more information.
The parsed metadata is used to make a counts image that is an idealized image containing the number of photons each WFI pixel would collect over an observation. It includes no systematic effects or noise beyond Poisson noise from the astronomical scene and backgrounds. Actual WFI observations are more complicated than just noisy versions of this idealized image, however, for several reasons:
WFI pixels have a very uncertain pedestal.
WFI pixels are sampled “up the ramp” during an observation, so a number of reads contribute to the final estimate for the rate of photons entering each pixel.
WFI reads are averaged on the telescope into resultants; ground images see only resultants.
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 here.
L2 images are currently constructed from count images at present simply by adding the following effects: reciprocity failure, nonlinearity, interpixel capacitance, and read noise. Then the gain is divided out and the image is quantized. The L2 images also do not yet remove the mean dark count rate. Most critically, this process ignores the details of ramp fitting and how that translates into uncertainties on each pixel, but at least the pieces are in place to include that sort of effect. The L2 images should not presently be taken seriously.
romanisim.image Module¶
Roman WFI simulator tool.
Based on galsim’s implementation of Roman image simulation. Uses galsim Roman modules for most of the real work.
Functions¶
|
Wrap a galsim simulated image with ASDF/roman_datamodel metadata. |
|
Simulate an image in a filter given a total count image. |
|
This routine kicks the tires on everything in this module. |
|
Simulate a sequence of observations on a field in different bandpasses. |
|
Simulate total counts in a single SCA. |
Making L1 images¶
An L1 (level 1) image is a “raw” image received from the detectors. The actual measurements made on the spacecraft consist of a number of non-destructive reads of the pixels of the H4RG detectors. These reads have independent read noise but because the pixels count the total number of photons having entered each pixel, the Poisson noise in different reads of the same pixel is correlated.
Because the telescope has limited bandwidth, every read is not transferred to ground stations. Instead, reads are averaged into “resultants” according to a specification called a MultiAccum table, and these resultants are transferred, archived, and analyzed. These resultants make up an L1 image, which romanisim simulates.
L1 images are created using an idealized counts image described here, which contains the number of photons each pixel of the detector would receive absent any instrumental systematics. To transform this into an L1 image, these counts must be apportioned into reads and averaged into resultants, and instrumental effects must be added.
This process proceeds by simulating each read, drawing the appropriate number of photons from the total number of photons for each read following a binomial distribution. These photons are added to a running sum that is then averaged into a resultant according to the MultiAccum table specification. This process requires drawing random numbers from the binomial distribution for every read of every pixel, and so can take on the order of a minute, but it allows detailed simulation of the statistics of the noise in each resultant together with their correlations. It also makes it straightforward to add various instrumental effects into the simulation accurately, since these usually apply to individual reads rather than to resultants (e.g., cosmic rays affect individual reads, and their affect on a resultant depends on the read in the resultant to which they apply).
After apportioning counts to resultants, systematic effects are added to the resultants. Presently only read noise is added. The read noise is averaged down like \(1/\sqrt{N}\), where \(N\) is the number of reads contributing to the resultant.
romanisim.l1 Module¶
Convert images into L1 images, with ramps.
We imagine starting with an image that gives the total number of counts from all Poisson processes (at least: sky, sources, dark current). We then need to redistribute these counts over the resultants of an L1 image.
The easiest thing to do, and probably where I should start, is to sample the image read-by-read with a binomial draw from the total counts weighted by the chance that each count landed in this particular time window. Then gather those for each resultant and average, as done on the spacecraft.
It’s tempting to go straight to making the appropriate resultants. Following Casertano (2022?), the variance in each resultant is:
where f is the count rate, N is the number of reads in the resultant, and \(\tau\) is the ‘variance-based resultant time’
where the t_k is the time of the kth read in the resultant.
For uniformly spaced reads,
where t_0 is the time of the first read in the resultant and d is the spacing of the reads.
So that gives the variance from Poisson counts in resultant. But how should we draw random numbers to get that variance and the right mean? I can separately control the mean and variance by scaling the Poisson distribution, but I’m not sure that’s doing the right thing with the higher order moments.
It probably isn’t that expensive to just work out all of the reads, individually, which will also allow more natural incorporation of cosmic rays down the road. So let’s take that approach instead for the moment.
How do we want to specify an L1 image? An L1 image is defined by a total count image and a list of lists \(t_{i, j}\), where \(t_{i, j}\) is the time at which the jth read in the ith resultant is made. We demand \(t_{i, j} > t_{k, l}\) whenever i > k or whenever i = k and j > l.
Things this doesn’t allow neatly:
jitter in telescope pointing: the rate image is the same for each read/resultant
non-linearity?
weird non-linear systematics in darks?
Possibly some systematics need to be applied to the individual reads, rather than to the final image. e.g., clearly nonlinearity? I need to think about when in the chain things like IPC, etc., come in. But it still seems correct to first generate the total number of counts that an ideal detector would measure from sources, and then apply these effects read-by-read as we build up the resultants—i.e., I expect the current framework will be able to handle this without major issue.
This approach is not super fast. For a high latitude set of resultants, generating all of the random numbers to determine the apportionment takes 43 s on the machine I’m currently using; this will scale linearly with the number of reads. That’s longer than the actual image production for the dummy scene I’m using (though only ~2x longer).
I don’t have a good way to speed this up. Explicitly doing the Poisson noise on each read from a rate image (rather than the apportionment of an image that already has Poisson noise) is 2x slower—generating billions of random numbers just takes a little while.
Plausibly I could figure out how to draw numbers directly from what a resultant is rather than looking at each read individually. That would likely bring a ~10x speed-up. The read noise there is easy. The poisson noise is a sum of scaled Poisson variables:
where \(c_i\) is a Poisson-distributed variable. The sum of Poisson-distributed variables is Poisson-distributed, but I wasn’t immediately able to find anything about the sum of scaled Poisson-distributed variables. The result is clearly not Poisson-distributed, but maybe there’s some special way to sample from that directly.
If we did sample from that directly, we’d still also need to get the sum of the counts in the reads comprising the resultant. So I think you’d need a separate draw for that, conditional on the number you got for the resultant. Or, reversing that, you might want to draw the total number of counts first, e.g., via the binomial distribution, and then you’d want to draw a number for what the average number of counts was among the reads comprising the resultant, conditional on the total number of counts. Then
is some kind of statistic of the multinomial distribution. That sounds a little more tractable?
We want to draw from \(\sum (N-i) c_i\). I think the probabilities are always \(1/N\), with the possible small but important exception of ‘skipped’ or ‘dropped’ reads, in which case the first read would be more like \(2/(N+1)\) and all the others \(1/(N+1)\). If the probabilities were always 1/N, this looks vaguely like it could have a nice analytic solution. Otherwise, I don’t immediately see a route forward. So I am not going to pursue this avenue further.
Functions¶
|
Adds read noise to resultants. |
|
Apportion counts to resultants given read times. |
|
Get the times of each read going into resultants for a MA table. |
|
Package and optionally write out an L1 frame. |
|
Make an L1 image from a counts image. |
|
Convert a set of times tij to corresponding probabilities for sampling. |
|
Verify that a set of times tij for a valid resultant. |
Catalogs¶
The simulator takes catalogs describing objects in a scene and generates images of that scene. These catalogs have the following form:
ra dec type n half_light_radius pa ba F087
float64 float64 str3 float64 float64 float64 float64 float64
------- ------- ---- ------- ----------------- ------- ------- --------
269.9 66.0 SER 1.6 0.6 165.6 0.9 1.80e-09
270.1 66.0 SER 3.6 0.4 71.5 0.7 3.35e-09
269.8 66.0 PSF -1.0 0.0 0.0 1.0 2.97e-10
269.9 66.0 SER 2.5 0.8 308.8 0.7 1.50e-09
269.8 65.9 SER 3.9 0.9 210.0 0.9 3.28e-10
270.1 66.0 SER 4.0 1.1 225.1 1.0 1.61e-09
269.9 65.9 SER 1.5 0.3 271.8 0.6 1.13e-09
269.9 65.9 SER 2.9 2.3 27.6 1.0 3.28e-09
269.9 66.0 SER 1.1 0.3 4.3 1.0 9.99e-10
The following fields must be specified for each source:
ra: the right ascension of the source
dec: the declination of the source
type: PSF or SER; whether the source is a point source or Sersic galaxy
n: the Sersic index of the source. This value is ignored if the source is a point source.
half_light_radius: the half light radius of the source in arcseconds. This value is ignored if the source is a point source.
pa: the position angle of the source, in degrees east of north. This value is ignored if the source is a point source.
ba: the major-to-minor axis ratio. This value is ignored if the source is a point source.
Following these required fields is a series of columns giving the fluxes of the the sources in “maggies”; the AB magnitude of the source is given by \(-2.5*\log_{10}(\mathrm{flux})\). In order to simulate a scene in a given bandpass, a column with the name of that bandpass must be present giving the total fluxes of the sources. Many flux columns may be present, and other columns may also be present but will be ignored.
The simulator then renders these images in the scene and produces the simulated L1 or L2 images.
romanisim.catalog Module¶
Catalog generation and reading routines.
This module provides basic routines to allow romanisim to render scenes based on catalogs of sources in those scenes.
Functions¶
|
Make a dummy catalog for testing purposes. |
|
Make a dummy table catalog. |
|
Read a catalog into a list of CatalogObjects. |
|
Read a astropy Table into a list of CatalogObjects. |
Classes¶
|
Simple class to hold galsim positions and profiles of objects. |
Class Inheritance Diagram¶
APT file support¶
The simulator possesses rudimentary support for simulating images from APT files. In order to simulate a scene, romanisim needs to know what’s in the scene, as specified by a catalog. It also needs to know where the telescope is pointed, the roll angle of the telescope, the date of the observation, and the bandpass. Finally, it needs to know what the MultiAccum table of the observation is—roughly, how long the exposure is and how the reads of the detector should be averaged into resultants.
Much of this information is available in an APT file. A rudimentary APT file reader can pull out the right ascension and declinations of observations, as well as the filters requested. However, support for roll angles is not yet included. APT files do not include the dates of observation, so this likewise is not included. APT files naturally do not contain catalogs of sources in the field, so some provision must be made for adding this information.
This module is not yet fully baked.
romanisim.apt Module¶
Very simple APT reader.
Converts an APT file into a list of (ra, dec, angle, filter, date, exposure time) needed for generating observations. This is adequate for reading in a few of the example Roman APTs but only supports a tiny fraction of what an APT file seems able to do.
Functions¶
|
Read an APT file, returning a list of observations. |
Classes¶
|
An observation of a target. |
|
A target for observation. |
Class Inheritance Diagram¶
Bandpasses¶
The simulator can render scenes in a number of different bandpasses. The choice of bandpass affects the point spread function used, the sky backgrounds, the fluxes of sources, and the reference files requested.
At present, romanisim simply passes the choice of bandpass to other packages—to webbpsf for PSF modeling, to galsim.roman for sky background estimation, to CRDS for reference file selection, or to the catalog for the selection of appropriate fluxes. However, because catalog fluxes are specified in “maggies” (i.e., in linear fluxes on the AB scale), the simulator needs to know how to convert between a maggie and the number of photons Roman receives from a source. Accordingly, the simulator knows about the AB zero points of the Roman filters, as derived from https://roman.gsfc.nasa.gov/science/WFI_technical.html .
romanisim.bandpass Module¶
Roman bandpass routines
The primary purpose of this module is to provide the number of counts per second expected for sources observed by Roman given a source with the nominal flat AB spectrum of 3631 Jy. The ultimate source of this information is https://roman.gsfc.nasa.gov/science/WFI_technical.html .
Functions¶
|
Compute the AB zero point fluxes for each filter. |
|
|
|
Read an effective area file from Roman. |
Point Spread Function Modeling¶
The simulator has two mechanisms for point source modeling. The first uses the galsim implementation of the Roman point spread function; for more information, see the galsim Roman documentation. The second uses the webbpsf package to make a model of the Roman PSF.
In the current implementation, the simulator uses a spatially constant, achromatic bandpass for each filter when using webbpsf. That is, the PSF is not modeled as varying across each SCA, nor does the PSF vary depending on the spectrum of the source being rendered. However, it seems straightforward to implement either of these modes in the context of galsim, albeit at some computational expense.
When using the galsim PSF, galsim’s “photon shooting” mode is used for efficient rendering of chromatic sources. When using webbpsf, FFTs are used to to do the convolution of the intrinsic source profile with the PSF and pixel grid of the instrument.
romanisim.psf Module¶
Roman PSF interface for galsim.
galsim.roman has an implementation of Roman’s PSF based on the aperture and some estimates for the wavefront errors over the aperture as described by amplitudes of various Zernicke modes. This seems like a very good approach, but we want to add here a mode using the official PSFs coming out of webbpsf, which takes a very similar overall approach.
galsim’s InterpolatedImage class makes this straightforward. Future work should consider the following:
how do we want to deal with PSF variation over the field? Can we be more efficient than making a new InterpolatedImage at each location based on some e.g. quadratic approximation to the PSF’s variation with location in an SCA?
how do we want to deal with the dependence of the PSF on the source SED? It’s possible I can just subclass ChromaticObject and implement evaluateAtWavelength, possibly also stealing the _shoot code from ChromaticOpticalPSF?
Functions¶
|
Make a PSF profile for Roman. |
World Coordinate Systems & Distortion¶
The simulator has two options for modeling distortion and world coordinate systems. The first is to use the routines in the galsim.roman package; see GalSim’s documentation for more information. The second is to use distortion reference files from the Calibration References Data System (CRDS).
The latter system works by grabbing reference distortion maps for the appropriate detector and filter combinations from the Roman CRDS server. These distortion maps are then wrapped into a galsim WCS object and fed to galsim’s rendering pipeline.
romanisim.wcs Module¶
Roman WCS interface for galsim.
galsim.roman has an implementation of Roman’s WCS based on some SIP coefficients for each SCA. This is presumably plenty good, but here we take the alternative approach of using the distortion functions provided in CRDS. These naturally are handled by the gWCS library, but galsim only naturally supports astropy WCSes via the ~legacy interface. So this module primarily makes the wrapper that interfaces gWCS and galsim.CelestialWCS together.
This presently gives rather different world coordinates given a specific telescope boresight. Partially this is not doing the same roll_ref determination that galsim.roman does, which could be fixed. But additionally the center of the SCA looks to be in a different place relative to the boresight for galsim.roman than for what I get from CRDS. This bears more investigation.
Functions¶
|
Get a WCS object for a given sca or set of CRDS parameters. |
|
Create a gWCS from a target position, a roll, and a distortion map. |
Classes¶
|
This WCS uses gWCS to implent a galsim CelestialWCS. |
Class Inheritance Diagram¶
Parameters¶
The parameters module contains a few useful constants describing the Roman telescope. At present these consist only of:
the default read noise when CRDS is not used
the number of border pixels
the specification of the read indices going into resultants for a fiducial L1 image
the read time.
romanisim.parameters Module¶
Parameters class storing a few useful constants for Roman simulations.
Utilities¶
The simulator utility module consists of miscellaneous utility routines intended to be of broad use. Present examples include:
turning astropy coordinates into galsim coordinates and vice-versa
making an RGB image from image slices
generating points at random in a region on the sky
flattening & unflattening nested dictionaries.
romanisim.util Module¶
Miscellaneous utility routines.
Functions¶
|
Turn a SkyCoord into a CelestialCoord. |
Convert a set of nested dictionaries into a flattened dictionary. |
|
|
Choose locations at random within radius of coord. |
|
|
|
Turn a CelestialCoord into a SkyCoord. |
Convert a flattened dictionary into a set of nested dictionaries. |