diff --git a/docs/index.rst b/docs/index.rst index 9209302c..e82074b5 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -15,6 +15,7 @@ Contents romanisim/image romanisim/l1 romanisim/l2 + romanisim/refs romanisim/catalog romanisim/apt diff --git a/docs/romanisim/refs.rst b/docs/romanisim/refs.rst new file mode 100644 index 00000000..fa2558b9 --- /dev/null +++ b/docs/romanisim/refs.rst @@ -0,0 +1,32 @@ +Reference files +=============== + +romanisim uses reference files from `CRDS `_ 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. diff --git a/romanisim/catalog.py b/romanisim/catalog.py index de5bdab7..54ffb668 100644 --- a/romanisim/catalog.py +++ b/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 diff --git a/romanisim/image.py b/romanisim/image.py index e8a2b165..b01c8f71 100644 --- a/romanisim/image.py +++ b/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,7 +248,7 @@ 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), @@ -224,7 +256,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 * 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 diff --git a/romanisim/l1.py b/romanisim/l1.py index 2e7c5b3f..ea56c455 100644 --- a/romanisim/l1.py +++ b/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 diff --git a/romanisim/ramp.py b/romanisim/ramp.py index 798eee9f..9ddbd77a 100644 --- a/romanisim/ramp.py +++ b/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. diff --git a/romanisim/util.py b/romanisim/util.py index fddc2c6d..3be53337 100644 --- a/romanisim/util.py +++ b/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):