From e6f942f94cae83b0c4c94e900790b206379be6b0 Mon Sep 17 00:00:00 2001 From: Eddie Schlafly Date: Mon, 15 Apr 2024 11:48:45 -0400 Subject: [PATCH] Allow SCA -1 to indicate that multiple images should be made. (#109) * Allow specifying SCA "negative one" to indicate that romanisim-make-image should make all 18 SCAs sequentially. --- .github/workflows/ci.yml | 3 +- .github/workflows/ci_cron.yml | 1 - pyproject.toml | 4 +- romanisim/image.py | 12 ++-- romanisim/psf.py | 3 +- romanisim/ris_make_utils.py | 16 ++++- romanisim/wcs.py | 19 ++++++ scripts/romancal_make_regtest_l1s.sh | 1 + scripts/romanisim-make-image | 92 +++++++++++++++++----------- 9 files changed, 99 insertions(+), 52 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 71e3c57f..297bc705 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -71,8 +71,7 @@ jobs: cache-path: ${{ needs.data.outputs.path }} cache-key: data-${{ needs.data.outputs.hash }} envs: | - - linux: py39-oldestdeps-cov-xdist - - linux: py39-xdist + - linux: py310-oldestdeps-cov-xdist - linux: py310-xdist - linux: py3-xdist - macos: py3-xdist diff --git a/.github/workflows/ci_cron.yml b/.github/workflows/ci_cron.yml index 4e594013..594b1a4d 100644 --- a/.github/workflows/ci_cron.yml +++ b/.github/workflows/ci_cron.yml @@ -60,6 +60,5 @@ jobs: cache-path: ${{ needs.data.outputs.path }} cache-key: data-${{ needs.data.outputs.hash }} envs: | - - macos: py39-xdist - macos: py310-xdist - linux: py3-devdeps-xdist diff --git a/pyproject.toml b/pyproject.toml index ac2933a5..3b50c905 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -2,7 +2,7 @@ name = 'romanisim' description = 'Nancy Grace Roman Space Telescope WFI Simulator' readme = 'README.md' -requires-python = '>=3.9' +requires-python = '>=3.10' license = { file = 'LICENSE' } authors = [{ name = 'STScI', email = 'help@stsci.edu' }] classifiers = [ @@ -13,7 +13,7 @@ classifiers = [ 'Operating System :: POSIX', 'Programming Language :: Python :: 3', 'Programming Language :: Python :: 3 :: Only', - 'Programming Language :: Python :: 3.9', + 'Programming Language :: Python :: 3.10', ] dependencies = [ 'asdf >=2.15.1', diff --git a/romanisim/image.py b/romanisim/image.py index 30063f44..43979b78 100644 --- a/romanisim/image.py +++ b/romanisim/image.py @@ -458,7 +458,9 @@ def simulate_counts_generic(image, exptime, objlist=None, psf=None, if not np.all(flat == 1): image.quantize() - image.array[:, :] = np.random.binomial( + rng_numpy_seed = rng.raw() + rng_numpy = np.random.default_rng(rng_numpy_seed) + image.array[:, :] = rng_numpy.binomial( image.array.astype('i4'), flat / maxflat) if dark is not None: @@ -800,6 +802,7 @@ def simulate(metadata, objlist, extras["simulate_reffiles"][key] = value extras['simcatobj'] = simcatobj + extras['wcs'] = wcs.convert_wcs_to_gwcs(counts.wcs) log.info('Simulation complete.') return im, extras @@ -878,12 +881,7 @@ def make_asdf(slope, slopevar_rn, slopevar_poisson, metadata=None, out['meta'].update(metadata) if imwcs is not None: # add a WCS - if isinstance(imwcs, wcs.GWCS): - out['meta'].update(wcs=imwcs.wcs) - else: - # make a gwcs WCS from a galsim.roman WCS - out['meta'].update( - wcs=wcs.wcs_from_fits_header(imwcs.header.header)) + out['meta'].update(wcs=wcs.convert_wcs_to_gwcs(imwcs)) out['data'] = slope out['dq'] = np.zeros(slope.shape, dtype='u4') diff --git a/romanisim/psf.py b/romanisim/psf.py index 28b84468..40294500 100644 --- a/romanisim/psf.py +++ b/romanisim/psf.py @@ -152,9 +152,10 @@ def make_psf(sca, filter_name, wcs=None, webbpsf=True, pix=None, pix=pix, chromatic=chromatic, **kw) elif pix is not None: raise ValueError('cannot set both pix and variable') - buf = 40 + buf = 49 # WebbPSF complains if we get too close to (0, 0) for some reason. # For other corners one can go to within a fraction of a pixel. + # if we go larger than 49 we have to change some of the tests, which use a 100x100 image. corners = dict( ll=[buf, buf], lr=[roman.n_pix - buf, buf], ul=[buf, roman.n_pix - buf], ur=[roman.n_pix - buf, roman.n_pix - buf]) diff --git a/romanisim/ris_make_utils.py b/romanisim/ris_make_utils.py index f9e88f95..a28e084e 100755 --- a/romanisim/ris_make_utils.py +++ b/romanisim/ris_make_utils.py @@ -4,6 +4,7 @@ from copy import deepcopy import os import re +import numpy as np import asdf from astropy import table from astropy import time @@ -13,7 +14,7 @@ import roman_datamodels from roman_datamodels import stnode from romanisim import catalog, image, wcs -from romanisim import parameters +from romanisim import parameters, log def merge_nested_dicts(dict1, dict2): @@ -148,6 +149,17 @@ def create_catalog(metadata=None, catalog_name=None, bandpasses=['F087'], coord, bandpasses=bandpasses, nobj=nobj, rng=rng) else: cat = table.Table.read(catalog_name, comment="#", delimiter=" ") + bandpass = [f for f in cat.dtype.names if f[0] == 'F'] + bad = np.zeros(len(cat), dtype='bool') + for b in bandpass: + bad |= ~np.isfinite(cat[b]) + if hasattr(cat[b], 'mask'): + bad |= cat[b].mask + cat = cat[~bad] + nbad = np.sum(bad) + if nbad > 0: + log.info(f'Removing {nbad} catalog entries with non-finite or ' + 'masked fluxes.') return cat @@ -251,7 +263,7 @@ def simulate_image_file(args, metadata, cat, rng=None, persist=None): args.pretend_spectral.upper()) drop_extra_dq = getattr(args, 'drop_extra_dq', False) - if drop_extra_dq: + if drop_extra_dq and ('dq' in romanisimdict): romanisimdict.pop('dq') # Write file diff --git a/romanisim/wcs.py b/romanisim/wcs.py index 2973f469..bd5b5e99 100644 --- a/romanisim/wcs.py +++ b/romanisim/wcs.py @@ -384,3 +384,22 @@ def coeffs_to_poly(mat, degree): gw.bounding_box = ((-0.5, nx - 0.5), (-0.5, ny - 0.5)) return gw + + +def convert_wcs_to_gwcs(wcs): + """Convert a GalSim WCS object into a GWCS object. + + Parameters + ---------- + wcs : gwcs.wcs.WCS or wcs.GWCS + input WCS to convert + + Returns + ------- + wcs.GWCS corresponding to wcs. + """ + if isinstance(wcs, GWCS): + return wcs.wcs + else: + # make a gwcs WCS from a galsim.roman WCS + return wcs_from_fits_header(wcs.header.header) diff --git a/scripts/romancal_make_regtest_l1s.sh b/scripts/romancal_make_regtest_l1s.sh index 6d35dec3..1bed276c 100644 --- a/scripts/romancal_make_regtest_l1s.sh +++ b/scripts/romancal_make_regtest_l1s.sh @@ -28,3 +28,4 @@ romanisim-make-image --radec 270.00 66.00 --level 1 --sca 1 --bandpass F158 --ca # truncated spectroscopic image romanisim-make-image --radec 270.00 66.00 --level 1 --sca 1 --bandpass F158 --catalog gaia-270-66-2027-06-01.ecsv --webbpsf --usecrds --ma_table_number 110 --date 2027-06-01T00:30:00 --rng_seed 7 --drop-extra-dq r0000201001001001001_01101_0004_WFI01_uncal.asdf --pretend-spectral GRISM & +wait diff --git a/scripts/romanisim-make-image b/scripts/romanisim-make-image index 63e68b64..c72ed21e 100755 --- a/scripts/romanisim-make-image +++ b/scripts/romanisim-make-image @@ -8,6 +8,60 @@ from astropy import units as u import galsim from romanisim import log, wcs, persistence, parameters from romanisim import ris_make_utils as ris +from copy import copy + + +def go(args): + if args.config is not None: + # Open and parse overrides file + with open(args.config, "r") as config_file: + config = yaml.safe_load(config_file) + 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 + + if args.sca == -1: + # simulate all 18 SCAs sequentially + origfilename = copy(args.filename) + origprevious = copy(args.previous) + for i in range(1, 19): + args.sca = i + args.filename = origfilename.format(f'WFI{args.sca:02d}') + if origprevious is not None: + args.previous = origprevious.format(f'WFI{args.sca:02d}') + go(args) + return + + 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, + frame='icrs') + wcs.fill_in_parameters(metadata, coord, boresight=args.boresight, pa_aper=args.roll) + + rng = galsim.UniformDeviate(args.rng_seed) + + # Create catalog + cat = ris.create_catalog(metadata=metadata, catalog_name=args.catalog, + bandpasses=[args.bandpass], rng=rng, + nobj=args.nobj, usecrds=args.usecrds) + + # Create persistence object + if args.previous is not None: + persist = persistence.Persistence.read(args.previous) + else: + persist = persistence.Persistence() + + # Simulate image and write to file + ris.simulate_image_file(args, metadata, cat, rng, persist) + if __name__ == '__main__': parser = argparse.ArgumentParser( @@ -61,40 +115,4 @@ if __name__ == '__main__': "packages like galsim's roman package or STIPS may better " "serve such purposes.") - if args.config is not None: - # Open and parse overrides file - with open(args.config, "r") as config_file: - config = yaml.safe_load(config_file) - 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, - truncate=args.truncate) - - if args.radec is not None: - coord = SkyCoord(ra=args.radec[0] * u.deg, dec=args.radec[1] * u.deg, - frame='icrs') - wcs.fill_in_parameters(metadata, coord, boresight=args.boresight, pa_aper=args.roll) - - rng = galsim.UniformDeviate(args.rng_seed) - - # Create catalog - cat = ris.create_catalog(metadata=metadata, catalog_name=args.catalog, - bandpasses=[args.bandpass], rng=rng, - nobj=args.nobj, usecrds=args.usecrds) - - # Create persistence object - if args.previous is not None: - persist = persistence.Persistence.read(args.previous) - else: - persist = persistence.Persistence() - - # Simulate image and write to file - ris.simulate_image_file(args, metadata, cat, rng, persist) + go(args)