From 30837260e82b38aff79db9b3163f2b25349131d8 Mon Sep 17 00:00:00 2001 From: Eddie Schlafly Date: Tue, 13 Feb 2024 13:20:55 -0500 Subject: [PATCH] Miscellaneous updates following RTB review (#99) * Updates following discussion with RTB. * Fix docs. --- docs/romanisim/overview.rst | 2 +- docs/romanisim/running.rst | 4 +- romanisim/catalog.py | 8 +- romanisim/image.py | 243 +++++++++++++++++++--------------- romanisim/l1.py | 13 +- romanisim/nonlinearity.py | 2 +- romanisim/parameters.py | 20 +-- romanisim/ramp.py | 3 - romanisim/ris_make_utils.py | 8 +- romanisim/tests/test_image.py | 19 +-- romanisim/tests/test_wcs.py | 1 + romanisim/wcs.py | 11 +- scripts/romanisim-make-image | 19 ++- scripts/romanisim-make-stack | 4 + 14 files changed, 205 insertions(+), 152 deletions(-) diff --git a/docs/romanisim/overview.rst b/docs/romanisim/overview.rst index 283c48e4..d053eee8 100644 --- a/docs/romanisim/overview.rst +++ b/docs/romanisim/overview.rst @@ -66,7 +66,7 @@ adds an additional top-level branch to the asdf tree with the name ├─ma_table_number (int): 1 ├─radec (NoneType): None ├─sca (int): 7 - ├─seed (NoneType): None + ├─rng_seed (NoneType): None ├─simcatobj (NDArrayType): shape=(496,), dtype=void96 ├─usecrds (bool): False └─webbpsf (bool): True diff --git a/docs/romanisim/running.rst b/docs/romanisim/running.rst index fd4302e2..42d59f81 100644 --- a/docs/romanisim/running.rst +++ b/docs/romanisim/running.rst @@ -39,7 +39,7 @@ this functionality:: microsecond (default: None) --level LEVEL 1 or 2, for L1 or L2 output (default: 2) --ma_table_number MA_TABLE_NUMBER - --seed SEED + --rng_seed SEED --nobj NOBJ --boresight radec specifies location of boresight, not center of WFI. (default: False) --previous PREVIOUS previous simulated file in chronological order used for persistence modeling. @@ -61,7 +61,7 @@ is available. The ``--webbpsf`` argument indicates that the `WebbPSF `_ package should be used to simulate the PSF; note that this presently disables chromatic PSF rendering. -The ``--seed`` argument specifies a seed to the random number +The ``--rng_seed`` argument specifies a seed to the random number generator, enabling reproducible results. The ``--nobj`` argument is only used when a catalog is not specified, diff --git a/romanisim/catalog.py b/romanisim/catalog.py index de49855e..77129c84 100644 --- a/romanisim/catalog.py +++ b/romanisim/catalog.py @@ -174,6 +174,9 @@ def make_galaxies(coord, power law index of magnitudes faintmag : float faintest AB magnitude for which to generate sources + Note this magnitude is in a "fiducial" band which is not observed. + Actual requested bandpasses are equal to this fiducial band plus + 1 mag of Gaussian noise. hlr_at_faintmag : float typical half light radius at faintmag (arcsec) bandpasses : list[str] @@ -269,6 +272,9 @@ def make_stars(coord, implies 3/5. faintmag : float faintest AB magnitude for which to generate sources + Note this magnitude is in a "fiducial" band which is not observed. + Actual requested bandpasses are equal to this fiducial band plus + 1 mag of Gaussian noise. truncation_radius : float truncation radius of cluster if not None; otherwise ignored. bandpasses : list[str] @@ -373,7 +379,7 @@ def table_to_catalog(table, bandpasses): obj = galsim.Sersic(table['n'][i], table['half_light_radius'][i]) obj = obj.shear( q=table['ba'][i], - beta=(table['pa'][i] + np.pi / 2) * galsim.radians) + beta=(np.radians(table['pa'][i]) + np.pi / 2) * galsim.radians) else: raise ValueError('Catalog types must be either PSF or SER.') out.append(CatalogObject(pos, obj, fluxes)) diff --git a/romanisim/image.py b/romanisim/image.py index d2f82f50..c5ef2636 100644 --- a/romanisim/image.py +++ b/romanisim/image.py @@ -65,7 +65,7 @@ def make_l2(resultants, read_pattern, read_noise=None, gain=None, flat=None, - linearity=None, dark=None, dq=None): + linearity=None, darkrate=None, dq=None): """ Simulate an image in a filter given resultants. @@ -83,8 +83,8 @@ def make_l2(resultants, read_pattern, read_noise=None, gain=None, flat=None, flat field to use linearity : romanisim.nonlinearity.NL object or None non-linearity correction to use. - dark : np.ndarray[nresultants, nx, ny] (float) - dark image to subtract from ramps (DN) + darkrate : np.ndarray[nx, ny] (float) + dark rate image to subtract from ramps (electron / s) dq : np.ndarray[nresultants, nx, ny] (int) DQ image corresponding to resultants @@ -99,19 +99,14 @@ def make_l2(resultants, read_pattern, read_noise=None, gain=None, flat=None, """ if read_noise is None: - read_noise = parameters.read_noise + read_noise = parameters.reference_data['readnoise'] if gain is None: - gain = parameters.gain + gain = parameters.reference_data['gain'] if linearity is not None: resultants = linearity.apply(resultants) - # no error propagation - - if dark is not None: - resultants = resultants - dark - log.info('Fitting ramps.') # commented out code below is inverse-covariance ramp fitting @@ -130,12 +125,11 @@ def make_l2(resultants, read_pattern, read_noise=None, gain=None, flat=None, ramppar, rampvar = ramp.fit_ramps_casertano( resultants * gain, dq & parameters.dq_do_not_use, 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 ' - 'a separate dark rate term to fit ramps.') # could iterate if we wanted to improve the flux estimates + if darkrate is not None: + ramppar[..., 1] -= darkrate + # The ramp fitter is not presently unit-aware; fix up the units by hand. # To do this right the ramp fitter should be made unit aware. # It takes a bit of work to get this right because we use the fact @@ -143,7 +137,7 @@ def make_l2(resultants, read_pattern, read_noise=None, gain=None, flat=None, # which isn't true as soon as things start having units and requires # special handling. And we use read_time without units a lot throughout # the code base. - slopes = ramppar[..., 1] / u.s + slopes = ramppar[..., 1] * u.electron / u.s readvar = rampvar[..., 0, 1, 1] * (u.electron / u.s)**2 poissonvar = rampvar[..., 1, 1, 1] * (u.electron / u.s)**2 @@ -556,6 +550,113 @@ def simulate_counts(metadata, objlist, return image, simcatobj +def gather_reference_data(image_mod, usecrds=False): + """Gather reference data corresponding to metadata. + + This function pulls files from parameters.reference_data and/or + CRDS to fill out the various reference files needed to perform + the simulation. If CRDS is set, values in parameters.reference_data + are used instead of CRDS files when the reference_data are None. If + all CRDS files should be used, parameters.reference_data must contain + only Nones. + + This functionality is intended to allow users to specify different + levels via a configuration file and not have them be overwritten + by the CRDS defaults, but it's not terribly clean. + + The input metadata is updated with CRDS software versions if CRDS + is used. + + Returns + ------- + dictionary containing the following keys: + read_noise + darkrate + gain + inv_linearity + linearity + saturation + reffiles + These have the reference images or constant values for the various + reference parameters. + """ + + reffiles = {k: v for k, v in parameters.reference_data.items()} + + out = dict(**reffiles) + if usecrds: + import crds + refsneeded = [k for (k, v) in reffiles.items() if v is None] + flatneeded = 'flat' in refsneeded + if flatneeded: + refsneeded.remove('flat') + if len(refsneeded) > 0: + reffiles.update(crds.getreferences( + image_mod.get_crds_parameters(), reftypes=refsneeded, + observatory='roman')) + if flatneeded: + try: + flatfile = crds.getreferences( + image_mod.get_crds_parameters(), + reftypes=['flat'], observatory='roman')['flat'] + + flat_model = roman_datamodels.datamodels.FlatRefModel(flatfile) + flat = flat_model.data[...].copy() + except crds.core.exceptions.CrdsLookupError: + log.warning('Could not find flat; using 1') + flat = 1 + out['flat'] = flat + image_mod.meta.ref_file.crds.sw_version = crds.__version__ + image_mod.meta.ref_file.crds.context_used = crds.get_context_name( + observatory=image_mod.crds_observatory) + + # reffiles has all of the reference files / values we know about + + nborder = parameters.nborder + + # we now need to extract the relevant fields + if isinstance(reffiles['readnoise'], str): + model = roman_datamodels.datamodels.ReadnoiseRefModel( + reffiles['readnoise']) + out['readnoise'] = model.data[nborder:-nborder, nborder:-nborder].copy() + + if isinstance(reffiles['gain'], str): + model = roman_datamodels.datamodels.GainRefModel(reffiles['gain']) + out['gain'] = model.data[nborder:-nborder, nborder:-nborder].copy() + + if isinstance(reffiles['dark'], str): + model = roman_datamodels.datamodels.DarkRefModel(reffiles['dark']) + out['dark'] = model.dark_slope[nborder:-nborder, nborder:-nborder].copy() + out['dark'] *= out['gain'] + if isinstance(out['dark'], u.Quantity): + out['dark'] = out['dark'].to(u.electron / u.s).value + + if isinstance(reffiles['inverselinearity'], str): + ilin_model = roman_datamodels.datamodels.InverselinearityRefModel( + reffiles['inverselinearity']) + out['inverselinearity'] = nonlinearity.NL( + ilin_model.coeffs[:, nborder:-nborder, nborder:-nborder].copy(), + ilin_model.dq[nborder:-nborder, nborder:-nborder].copy(), + gain=out['gain']) + + if isinstance(reffiles['linearity'], str): + lin_model = roman_datamodels.datamodels.LinearityRefModel( + reffiles['linearity']) + out['linearity'] = nonlinearity.NL( + lin_model.coeffs[:, nborder:-nborder, nborder:-nborder].copy(), + lin_model.dq[nborder:-nborder, nborder:-nborder].copy(), + gain=out['gain']) + + if isinstance(reffiles['saturation'], str): + saturation = roman_datamodels.datamodels.SaturationRefModel( + reffiles['saturation']) + saturation = saturation.data[nborder:-nborder, nborder:-nborder].copy() + out['saturation'] = saturation + + out['reffiles'] = reffiles + return out + + def simulate(metadata, objlist, usecrds=True, webbpsf=True, level=2, crparam=dict(), persistence=None, seed=None, rng=None, **kwargs @@ -603,17 +704,23 @@ def simulate(metadata, objlist, also include information on persistence and cosmic ray hits. """ + if not usecrds: + log.warning('--usecrds is not set. romanisim will not use reference ' + 'files from CRDS. The WCS may be incorrect and up-to-date ' + 'calibration information will not be used.') + meta = maker_utils.mk_common_meta() meta["photometry"] = maker_utils.mk_photometry() + meta['wcs'] = None for key in parameters.default_parameters_dictionary.keys(): meta[key].update(parameters.default_parameters_dictionary[key]) - util.add_more_metadata(meta) - for key in metadata.keys(): meta[key].update(metadata[key]) + util.add_more_metadata(meta) + # Create Image model to track validation image_node = maker_utils.mk_level2_image() image_node['meta'] = meta @@ -623,86 +730,16 @@ def simulate(metadata, objlist, filter_name = image_mod.meta.instrument.optical_element 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 - # data in numpy arrays. - # should query CRDS for any reference files not specified on the command - # line. - reffiles = {} - if usecrds: - import crds - refnames_lst = ['readnoise', 'dark', 'gain', 'inverselinearity', 'linearity', 'saturation'] - - for key, value in parameters.reference_data.items(): - if (key in refnames_lst) and (value is not None): - reffiles[key] = value - refnames_lst.remove(key) - - if refnames_lst: - reffiles.update(crds.getreferences( - image_mod.get_crds_parameters(), reftypes=refnames_lst, - observatory='roman')) - - read_noise_model = roman_datamodels.datamodels.ReadnoiseRefModel( - reffiles['readnoise']) - dark_model = roman_datamodels.datamodels.DarkRefModel( - reffiles['dark']) - gain_model = roman_datamodels.datamodels.GainRefModel( - reffiles['gain']) - inv_linearity_model = roman_datamodels.datamodels.InverselinearityRefModel( - reffiles['inverselinearity']) - linearity_model = roman_datamodels.datamodels.LinearityRefModel( - reffiles['linearity']) - saturation_model = roman_datamodels.datamodels.SaturationRefModel( - reffiles['saturation']) - try: - flatfile = crds.getreferences( - image_mod.get_crds_parameters(), - reftypes=['flat'], observatory='roman')['flat'] - - flat_model = roman_datamodels.datamodels.FlatRefModel(flatfile) - flat = flat_model.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. - nborder = parameters.nborder - read_noise = read_noise_model.data[nborder:-nborder, nborder:-nborder] - darkrate = dark_model.data[-1, nborder:-nborder, nborder:-nborder] / exptime_tau - dark = dark_model.data[:, nborder:-nborder, nborder:-nborder] - gain = gain_model.data[nborder:-nborder, nborder:-nborder] - inv_linearity = nonlinearity.NL( - inv_linearity_model.coeffs[:, nborder:-nborder, nborder:-nborder], - inv_linearity_model.dq[nborder:-nborder, nborder:-nborder], - gain=gain) - linearity = nonlinearity.NL( - linearity_model.coeffs[:, nborder:-nborder, nborder:-nborder], - linearity_model.dq[nborder:-nborder, nborder:-nborder], - gain=gain) - darkrate *= gain - image_mod.meta.ref_file.crds.sw_version = crds.__version__ - image_mod.meta.ref_file.crds.context_used = crds.get_context_name( - observatory=image_mod.crds_observatory - ) - saturation = saturation_model.data[nborder:-nborder, nborder:-nborder] - else: - read_noise = parameters.reference_data["readnoise"] \ - if parameters.reference_data["readnoise"] is not None else galsim.roman.read_noise - darkrate = parameters.reference_data["darkcurrent"] \ - if parameters.reference_data["darkcurrent"] is not None else galsim.roman.dark_current - dark = parameters.reference_data["dark"] - gain = parameters.reference_data["gain"] - flat = parameters.reference_data["flat"] \ - if parameters.reference_data["flat"] is not None else 1 - inv_linearity = parameters.reference_data["inverselinearity"] - linearity = parameters.reference_data["linearity"] - saturation = parameters.reference_data["saturation"] + refdata = gather_reference_data(image_mod, usecrds=usecrds) + read_noise = refdata['readnoise'] + darkrate = refdata['dark'] + gain = refdata['gain'] + inv_linearity = refdata['inverselinearity'] + linearity = refdata['linearity'] + saturation = refdata['saturation'] + reffiles = refdata['reffiles'] + flat = refdata['flat'] if rng is None and seed is None: seed = 43 @@ -735,7 +772,7 @@ def simulate(metadata, objlist, elif level == 2: slopeinfo = make_l2(l1, read_pattern, read_noise=read_noise, gain=gain, flat=flat, linearity=linearity, - dark=dark, dq=l1dq) + darkrate=darkrate, dq=l1dq) l2dq = np.bitwise_or.reduce(l1dq, axis=0) im, extras = make_asdf( *slopeinfo, metadata=image_mod.meta, persistence=persistence, @@ -792,10 +829,7 @@ def make_asdf(slope, slopevar_rn, slopevar_poisson, metadata=None, # cal['step'] is left empty for now, but in principle # could be filled out at some level # ephemeris contains a lot of angles that could be computed. - # exposure contains start_time, mid_time, end_time - # plus MJD equivalents - # plus TDB equivalents (?) - # plus ENG equivalents (?) + # exposure contains # ngroups, nframes, sca_number, gain_factor, integration_time, # elapsed_exposure_time, nints, integration_start, integration_end, # frame_divisor, groupgap, nsamples, sample_time, frame_time, @@ -825,13 +859,6 @@ def make_asdf(slope, slopevar_rn, slopevar_poisson, metadata=None, # leave border pixels blank for now. # amp33: leave blank for now. # data: image - # dq: I guess all 0? - # err: I guess I need to look up how this is defined. But - # the "right" thing is probably - # sqrt(noisless image + read_noise**2) - # i.e., Gaussian approximation of Poisson noise + read noise. - # this is just sqrt(var_poisson + var_rnoise + var_flat)? - # var_rnoise: okay # var_flat: currently zero if metadata is not None: out['meta'].update(metadata) diff --git a/romanisim/l1.py b/romanisim/l1.py index 3875f419..70090fd9 100644 --- a/romanisim/l1.py +++ b/romanisim/l1.py @@ -365,11 +365,14 @@ def add_read_noise_to_resultants(resultants, tij, read_noise=None, rng=None, else: rng = galsim.GaussianDeviate(rng) - noise = np.zeros(resultants.shape, dtype='f4') - rng.generate(noise) if read_noise is None: - read_noise = parameters.read_noise + read_noise = parameters.reference_data['readnoise'] + if read_noise is None: + log.warning('Not applying read noise due to weird reference data.') + return resultants + noise = np.zeros(resultants.shape, dtype='f4') + rng.generate(noise) noise = noise * read_noise / np.array( [len(x)**0.5 for x in tij]).reshape(-1, 1, 1) resultants += noise @@ -529,7 +532,7 @@ def make_l1(counts, read_pattern, resultants *= u.electron if gain is None: - gain = parameters.gain + gain = parameters.reference_data['gain'] if gain is not None and not isinstance(gain, u.Quantity): gain = gain * u.electron / u.DN log.warning('Making up units for gain.') @@ -554,7 +557,7 @@ def make_l1(counts, read_pattern, resultants += parameters.pedestal if saturation is None: - saturation = parameters.saturation + saturation = parameters.reference_data['saturation'] # this maybe should be better applied at read time? # it's not actually clear to me what the right thing to do diff --git a/romanisim/nonlinearity.py b/romanisim/nonlinearity.py index e6e4b2ee..ff603cc0 100644 --- a/romanisim/nonlinearity.py +++ b/romanisim/nonlinearity.py @@ -125,7 +125,7 @@ def __init__(self, coeffs, dq=None, gain=None): if dq is None: dq = np.zeros(coeffs.shape[1:], dtype='uint32') if gain is None: - gain = parameters.gain.to(u.electron / u.DN).value + gain = parameters.reference_data['gain'].to(u.electron / u.DN).value self.coeffs, self.dq = repair_coefficients(coeffs, dq) self.gain = gain diff --git a/romanisim/parameters.py b/romanisim/parameters.py index 2982d181..172a5e9d 100644 --- a/romanisim/parameters.py +++ b/romanisim/parameters.py @@ -41,25 +41,22 @@ 'v3_ref': 0, 'roll_ref': 0, }, + 'aperture': {'name': 'WFI_CEN', + 'position_angle': 0 + }, } reference_data = { - "dark": None, - "darkcurrent": None, + "dark": 0.01 * u.electron / u.s, "distortion": None, "flat": None, - "gain": None, + "gain": 2 * u.electron / u.DN, "inverselinearity": None, "linearity": None, - "readnoise": None, - "saturation": None + "readnoise": 5.0 * u.DN, + "saturation": 55000 * u.DN, } -# default read noise -read_noise = 5.0 * u.DN - -gain = 1 * u.electron / u.DN - nborder = 4 # number of border pixels used for reference pixels. read_time = 3.04 @@ -84,9 +81,6 @@ persistence = dict(A=0.017, x0=6.0e4, dx=5.0e4, alpha=0.045, gamma=1, half_well=50000, ignorerate=0.01) -# default saturation level in DN absent reference file information -saturation = 55000 * u.DN - # arbitrary constant to add to initial L1 image so that pixels aren't clipped at zero. pedestal = 100 * u.DN diff --git a/romanisim/ramp.py b/romanisim/ramp.py index 5ba7c840..caf94ed3 100644 --- a/romanisim/ramp.py +++ b/romanisim/ramp.py @@ -557,9 +557,6 @@ def fit_ramps_casertano(resultants, dq, read_noise, read_pattern): par = par[0] var = var[0] - if resultants_unit is not None: - par = par * resultants_unit - return par, var diff --git a/romanisim/ris_make_utils.py b/romanisim/ris_make_utils.py index e8cb1d31..0da9c7e6 100755 --- a/romanisim/ris_make_utils.py +++ b/romanisim/ris_make_utils.py @@ -38,7 +38,7 @@ def merge_nested_dicts(dict1, dict2): dict1[key] = value -def set_metadata(meta=None, date=None, bandpass='F087', sca=7, ma_table_number=1): +def set_metadata(meta=None, date=None, bandpass='F087', sca=7, ma_table_number=1, truncate=None): """ Set / Update metadata parameters @@ -75,6 +75,12 @@ def set_metadata(meta=None, date=None, bandpass='F087', sca=7, ma_table_number=1 # Observational metadata meta['instrument']['optical_element'] = bandpass meta['exposure']['ma_table_number'] = ma_table_number + meta['exposure']['read_pattern'] = parameters.read_pattern[ma_table_number] + if truncate is not None: + meta['exposure']['read_pattern'] = meta['exposure']['read_pattern'][:truncate] + meta['exposure']['truncated'] = True + else: + meta['exposure']['truncated'] = False return meta diff --git a/romanisim/tests/test_image.py b/romanisim/tests/test_image.py index 128453e1..8aff9e6a 100644 --- a/romanisim/tests/test_image.py +++ b/romanisim/tests/test_image.py @@ -75,37 +75,37 @@ def test_make_l2(): [9 + x for x in range(4)], [13 + x for x in range(4)]] slopes, readvar, poissonvar = image.make_l2( - resultants, read_pattern, gain=1, flat=1, dark=0) + resultants, read_pattern, gain=1, flat=1, darkrate=0) assert np.allclose(slopes, 0) resultants[:, :, :] = np.arange(4)[:, None, None] resultants *= u.DN gain = 1 * u.electron / u.DN slopes, readvar, poissonvar = image.make_l2( resultants, read_pattern, - gain=gain, flat=1, dark=0) + gain=gain, flat=1, darkrate=0) assert np.allclose(slopes, 1 / parameters.read_time / 4 * u.electron / u.s) assert np.all(np.array(slopes.shape) == np.array(readvar.shape)) assert np.all(np.array(slopes.shape) == np.array(poissonvar.shape)) assert np.all(readvar >= 0) assert np.all(poissonvar >= 0) slopes1, readvar1, poissonvar1 = image.make_l2( - resultants, read_pattern, read_noise=1, dark=0) + resultants, read_pattern, read_noise=1, darkrate=0) slopes2, readvar2, poissonvar2 = image.make_l2( - resultants, read_pattern, read_noise=2, dark=0) + resultants, read_pattern, read_noise=2, darkrate=0) assert np.all(readvar2 >= readvar1) # because we change the weights depending on the ratio of read & poisson # noise, we can't assume above that readvar2 = readvar1 * 4. # But it should be pretty darn close here. assert np.all(np.abs(readvar2 / (readvar1 * 4) - 1) < 0.1) slopes2, readvar2, poissonvar2 = image.make_l2( - resultants, read_pattern, read_noise=1, flat=0.5) + resultants, read_pattern, read_noise=1, flat=0.5, darkrate=0) assert np.allclose(slopes2, slopes1 / 0.5) assert np.allclose(readvar2, readvar1 / 0.5**2) assert np.allclose(poissonvar2, poissonvar1 / 0.5**2) slopes, readvar, poissonvar = image.make_l2( resultants, read_pattern, read_noise=1, gain=gain, flat=1, - dark=resultants) - assert np.allclose(slopes, 0) + darkrate=1 / parameters.read_time / 4) + assert np.allclose(slopes, 0, atol=1e-6) def set_up_image_rendering_things(): @@ -559,8 +559,9 @@ def test_simulate(): parameters.persistence['x0'], parameters.persistence['dx'], parameters.persistence['alpha'], parameters.persistence['gamma']) roughguess = roughguess * 140 # seconds of integration + gain = parameters.reference_data['gain'] assert np.abs( - np.log(np.mean(diff * parameters.gain).value / roughguess)) < 1 + np.log(np.mean(diff * gain).value / roughguess)) < 1 # within a factor of e log.info('DMS224: added persistence to an image.') @@ -723,7 +724,7 @@ def test_inject_source_into_image(): # Make new image of the combination newimage, readvar, poissonvar = image.make_l2( newramp * u.DN, read_pattern, - gain=1 * u.electron / u.DN, flat=1, dark=0) + gain=1 * u.electron / u.DN, flat=1, darkrate=0) # Create mask of PSF nonzero = (sourcecounts.array != 0) diff --git a/romanisim/tests/test_wcs.py b/romanisim/tests/test_wcs.py index caa5d5a5..96236f44 100644 --- a/romanisim/tests/test_wcs.py +++ b/romanisim/tests/test_wcs.py @@ -138,6 +138,7 @@ def test_wcs_crds_match(): for key in metadata.keys(): meta[key].update(metadata[key]) + meta['wcs'] = None image_node = maker_utils.mk_level2_image() image_node['meta'] = meta diff --git a/romanisim/wcs.py b/romanisim/wcs.py index 597caa95..9808fec3 100644 --- a/romanisim/wcs.py +++ b/romanisim/wcs.py @@ -44,8 +44,8 @@ def fill_in_parameters(parameters, coord, pa_aper=0, boresight=True): Parameters ---------- parameters : dict - CRDS parameters dictionary - keys like pointing.* and wcsinfo.* may be modified + Metadata dictionary + Dictionaries like pointing, aperture, and wcsinfo may be modified coord : astropy.coordinates.SkyCoord or galsim.CelestialCoord world coordinates at V2 / V3 ref (boresight or center of WFI CCDs) @@ -66,6 +66,8 @@ def fill_in_parameters(parameters, coord, pa_aper=0, boresight=True): if 'wcsinfo' not in parameters.keys(): parameters['wcsinfo'] = {} + if 'aperture' not in parameters.keys(): + parameters['aperture'] = {} parameters['wcsinfo']['ra_ref'] = ( parameters['pointing']['ra_v1']) @@ -81,10 +83,15 @@ def fill_in_parameters(parameters, coord, pa_aper=0, boresight=True): if boresight: parameters['wcsinfo']['v2_ref'] = 0 parameters['wcsinfo']['v3_ref'] = 0 + parameters['aperture']['name'] = 'BORESIGHT' else: from .parameters import v2v3_wficen parameters['wcsinfo']['v2_ref'] = v2v3_wficen[0] parameters['wcsinfo']['v3_ref'] = v2v3_wficen[1] + parameters['aperture']['name'] = 'WFI_CEN' + + parameters['aperture']['position_angle'] = pa_aper + def get_wcs(image, usecrds=True, distortion=None): diff --git a/scripts/romanisim-make-image b/scripts/romanisim-make-image index 50006c30..e4f2a97d 100755 --- a/scripts/romanisim-make-image +++ b/scripts/romanisim-make-image @@ -12,15 +12,15 @@ from romanisim import ris_make_utils as ris if __name__ == '__main__': parser = argparse.ArgumentParser( description='Make a demo image.', - epilog='EXAMPLE: %(prog)s output_image.fits', + epilog='EXAMPLE: %(prog)s output_image.asdf', formatter_class=argparse.ArgumentDefaultsHelpFormatter) - parser.add_argument('filename', type=str, help='output image (fits)') + parser.add_argument('filename', type=str, help='output image (asdf)') parser.add_argument('--bandpass', type=str, help='bandpass to simulate', default='F087') parser.add_argument('--boresight', action='store_true', default=False, help=('radec specifies location of boresight, not ' 'center of WFI.')) - parser.add_argument('--catalog', type=str, help='input catalog (csv)', + parser.add_argument('--catalog', type=str, help='input catalog (ecsv)', default=None) parser.add_argument('--config', type=str, help='input parameter override file (yaml)', default=None) @@ -43,6 +43,8 @@ if __name__ == '__main__': help='Use CRDS for distortion map') parser.add_argument('--webbpsf', action='store_true', help='Use webbpsf for PSF') + parser.add_argument('--truncate', type=int, default=None, help=( + 'If set, truncate the MA table at given number of resultants.')) args = parser.parse_args() @@ -61,10 +63,15 @@ if __name__ == '__main__': combo_dict = parameters.__dict__ ris.merge_nested_dicts(combo_dict, config) parameters.__dict__.update(combo_dict) + elif args.usecrds: + # don't use default values + for k in parameters.reference_data: + parameters.reference_data[k] = None - metadata = ris.set_metadata(date=args.date, bandpass=args.bandpass, - sca=args.sca, - ma_table_number=args.ma_table_number) + metadata = ris.set_metadata( + date=args.date, bandpass=args.bandpass, + sca=args.sca, ma_table_number=args.ma_table_number, + truncate=args.truncate) if args.radec is not None: coord = SkyCoord(ra=args.radec[0] * u.deg, dec=args.radec[1] * u.deg, diff --git a/scripts/romanisim-make-stack b/scripts/romanisim-make-stack index 66c833e2..582cc2cf 100755 --- a/scripts/romanisim-make-stack +++ b/scripts/romanisim-make-stack @@ -117,6 +117,10 @@ def main(): combo_dict = parameters.__dict__ ris.merge_nested_dicts(combo_dict, config) parameters.__dict__.update(combo_dict) + elif args.usecrds: + # don't use default values + for k in parameters.reference_data: + parameters.reference_data[k] = None # Set up metadata metadata = deepcopy(default_parameters_dictionary)