Skip to content

Commit

Permalink
Record simulation input arguments in simulation output.
Browse files Browse the repository at this point in the history
Record simulated objects' information in simulation output.
Exclude objects off frame more efficiently.
Fix exposure time scaling.
  • Loading branch information
schlafly committed Oct 18, 2022
1 parent 0f6391f commit 7a1fc56
Show file tree
Hide file tree
Showing 4 changed files with 89 additions and 21 deletions.
28 changes: 27 additions & 1 deletion docs/romanisim/overview.rst
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,33 @@ 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.
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.

The output of ``romanisim-make-image`` is an appropriate asdf file for
the requested level image, with the following addition. The script
adds an additional top-level branch to the asdf tree with the name
``romanisim``. Here's an example::

└─romanisim (dict)
├─bandpass (str): F087
├─catalog (NoneType): None
├─date (NoneType): None
├─filename (str): out.asdf
├─level (int): 1
├─ma_table_number (int): 1
├─radec (NoneType): None
├─sca (int): 7
├─seed (NoneType): None
├─simcatobj (NDArrayType): shape=(496,), dtype=void96
├─usecrds (bool): False
└─webbpsf (bool): True

These fields are simply the arguments to ``romanisim-make-image``,
plus an additional ``simcatobj`` field with contains the ``x``, ``y``,
and number of photons of each simulated source.

Features not included so far:

Expand Down
2 changes: 1 addition & 1 deletion romanisim/catalog.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ class CatalogObject:
"""Simple class to hold galsim positions and profiles of objects."""
sky_pos: galsim.CelestialCoord
profile: galsim.GSObject
flux : dict
flux: dict


def make_dummy_catalog(coord, radius=0.1, rng=None, seed=42, nobj=1000,
Expand Down
65 changes: 51 additions & 14 deletions romanisim/image.py
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,38 @@ def make_l2(resultants, ma_table, read_noise=None):
)


def in_bounds(objlist, wcs, imbd, margin):
"""Filter sources in an objlist to those landing on an image.
gWCS objects can be somewhat expensive to evaluate source-by-source, so
this transforms the object list back into an array to do this in a
vectorized way, and gives the results.
Parameters
----------
objlist : list[romanisim.catalog.CatalogObject]
list of objects
wcs : galsim.wcs.CelestialWCS
wcs of chip
imbd : galsim.Image.Bounds
bounds of image
margin : int
keep sources up to margin outside of bounds
Returns
-------
keep : np.ndarray (bool)
whether each source lands near the image (True) or not (False)
xx, yy : np.ndarray (float)
x, y locations of sources on chip
"""
coord = np.array([[o.sky_pos.ra.rad, o.sky_pos.dec.rad] for o in objlist])
xx, yy = wcs._xy(coord[:, 0], coord[:, 1])
keep = ((xx > imbd.xmin - margin) & (xx < imbd.xmax + margin) &
(yy > imbd.ymin - margin) & (yy < imbd.ymax + margin))
return keep, xx, yy


def simulate_counts(sca, targ_pos, date, objlist, filter_name,
exptime=None, rng=None, seed=None,
ignore_distant_sources=10, usecrds=True,
Expand Down Expand Up @@ -131,10 +163,12 @@ def simulate_counts(sca, targ_pos, date, objlist, filter_name,
Returns
-------
galsim.Image
image : galsim.Image
idealized image of scene as seen by Roman, giving total electron counts
from rate sources (astronomical objects; backgrounds; dark current) in
each pixel.
simcatobj : np.ndarray
catalog of simulated objects in image
"""
if exptime is None:
exptime = roman.exptime
Expand All @@ -145,6 +179,9 @@ def simulate_counts(sca, targ_pos, date, objlist, filter_name,
if rng is None:
rng = galsim.UniformDeviate(seed)

simcatobj = np.zeros(
len(objlist), dtype=[('x', 'f4'), ('y', 'f4'), ('photons', 'f4')])

galsim_filter_name = romanisim.bandpass.roman2galsim_bandpass[filter_name]
bandpass = roman.getBandpasses(AB_zeropoint=True)[galsim_filter_name]
imwcs = wcs.get_wcs(world_pos=targ_pos, date=date, sca=sca, usecrds=usecrds)
Expand All @@ -170,19 +207,15 @@ def simulate_counts(sca, targ_pos, date, objlist, filter_name,
nrender = 0
final = None
info = []
keep, xpos, ypos = in_bounds(objlist, imwcs, imbd, ignore_distant_sources)
for i, obj in enumerate(objlist):
# this is kind of slow. We need to do an initial vectorized cull before
# reaching this point.
t0 = time.time()
image_pos = imwcs.toImage(obj.sky_pos)
if ((image_pos.x < imbd.xmin - ignore_distant_sources) or
(image_pos.x > imbd.xmax + ignore_distant_sources) or
(image_pos.y < imbd.ymin - ignore_distant_sources) or
(image_pos.y > imbd.ymax + ignore_distant_sources)):
if not keep[i]:
# ignore source off edge. Could do better by looking at
# source size.
info.append(0)
continue
image_pos = galsim.PositionD(xpos[i], ypos[i])
final = galsim.Convolve(obj.profile * exptime, psf)
if chromatic:
stamp = final.drawImage(
Expand All @@ -191,7 +224,7 @@ def simulate_counts(sca, targ_pos, date, objlist, filter_name,
else:
if obj.flux is not None:
final = final.withFlux(
obj.flux[filter_name] * abflux)
obj.flux[filter_name] * abflux * exptime)
try:
stamp = final.drawImage(center=image_pos,
wcs=imwcs.local(image_pos))
Expand All @@ -201,10 +234,12 @@ def simulate_counts(sca, targ_pos, date, objlist, filter_name,
bounds = stamp.bounds & full_image.bounds
if bounds.area() == 0:
continue
simcatobj[i] = (xpos[i], ypos[i], np.sum(stamp.array))
full_image[bounds] += stamp[bounds]
nrender += 1
info.append(time.time() - t0)
log.info('Rendered %d sources...' % nrender)
simcatobj = simcatobj[keep]

poisson_noise = galsim.PoissonNoise(rng)
sky_image.addNoise(poisson_noise)
Expand All @@ -226,7 +261,7 @@ def simulate_counts(sca, targ_pos, date, objlist, filter_name,
if return_info:
full_image = (full_image, info)

return full_image
return full_image, simcatobj


def simulate(metadata, objlist,
Expand Down Expand Up @@ -262,9 +297,11 @@ def simulate(metadata, objlist,
Returns
-------
asdf structure with simulated image
image : roman_datamodels model
simulated image
simcatobj : np.ndarray
image positions and fluxes of simulated objects
"""

all_metadata = copy.deepcopy(parameters.default_parameters_dictionary)
flatmetadata = util.flatten_dictionary(metadata)
flatmetadata = {'roman.meta'+k if k.find('roman.meta') != 0 else k: v
Expand Down Expand Up @@ -309,7 +346,7 @@ def simulate(metadata, objlist,


log.info('Simulating filter {0}...'.format(filter_name))
counts = simulate_counts(
counts, simcatobj = simulate_counts(
sca, coord, date, objlist, filter_name, rng=rng,
exptime=exptime, usecrds=usecrds, darkrate=darkrate,
webbpsf=webbpsf)
Expand All @@ -322,7 +359,7 @@ def simulate(metadata, objlist,
elif level == 2:
slopeinfo = make_l2(l1, ma_table, read_noise=read_noise)
im = make_asdf(*slopeinfo, metadata=all_metadata)
return im
return im, simcatobj


def make_test_catalog_and_images(seed=12345, sca=7, filters=None, nobj=1000,
Expand Down
15 changes: 10 additions & 5 deletions scripts/romanisim-make-image
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,8 @@ from romanisim import log
if __name__ == '__main__':
parser = argparse.ArgumentParser(
description='Make a demo image.',
epilog='EXAMPLE: %(prog)s output_image.fits')
epilog='EXAMPLE: %(prog)s output_image.fits',
formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser.add_argument('filename', type=str, help='output image (fits)')
parser.add_argument('--catalog', type=str, help='input catalog (csv)',
default=None)
Expand Down Expand Up @@ -75,10 +76,14 @@ if __name__ == '__main__':
cat = Table.read(args.catalog)
objlist = catalog.table_to_catalog(cat, [args.bandpass])
rng = galsim.UniformDeviate(args.seed)
im = image.simulate(metadata, objlist, usecrds=args.usecrds,
webbpsf=args.webbpsf, level=args.level,
rng=rng)
im, simcatobj = image.simulate(
metadata, objlist, usecrds=args.usecrds,
webbpsf=args.webbpsf, level=args.level,
rng=rng)

romanisimdict = vars(args)
romanisimdict['simcatobj'] = simcatobj

af = asdf.AsdfFile()
af.tree = {'roman': im}
af.tree = {'roman': im, 'romanisim': romanisimdict}
af.write_to(args.filename)

0 comments on commit 7a1fc56

Please sign in to comment.