From 8dbcd8da4105234006aeb6dfcb5ee08a1281c32f Mon Sep 17 00:00:00 2001 From: PaulHuwe Date: Mon, 24 Jun 2024 13:25:44 -0400 Subject: [PATCH 1/9] Initial save. --- romanisim/image.py | 5 +- romanisim/l3.py | 396 +++++++++++++++++++++++++++++++++++-- romanisim/parameters.py | 9 + romanisim/tests/test_l3.py | 160 +++++++++++++++ romanisim/util.py | 84 ++++++++ romanisim/wcs.py | 16 +- 6 files changed, 649 insertions(+), 21 deletions(-) create mode 100644 romanisim/tests/test_l3.py diff --git a/romanisim/image.py b/romanisim/image.py index 601f442c..0727f3e5 100644 --- a/romanisim/image.py +++ b/romanisim/image.py @@ -29,10 +29,7 @@ from romanisim import log import roman_datamodels -try: - import roman_datamodels.maker_utils as maker_utils -except ImportError: - import roman_datamodels.testing.utils as maker_utils +import roman_datamodels.maker_utils as maker_utils # galsim fluxes are in photons / cm^2 / s diff --git a/romanisim/l3.py b/romanisim/l3.py index 0b47022c..db5a5af0 100644 --- a/romanisim/l3.py +++ b/romanisim/l3.py @@ -1,11 +1,29 @@ -"""Function for Level 3-like images. +"""Roman WFI simulator functions for Level 3 mosaics. +Based on galsim's implementation of Roman image simulation. Uses galsim Roman modules +for most of the real work. """ import numpy as np +import math +import astropy.time +from astropy import table import galsim -from . import log -from . import image, wcs, psf +from galsim import roman + +from . import parameters +from . import util +import romanisim.wcs +import romanisim.l1 +import romanisim.bandpass +import romanisim.psf +import romanisim.image +import romanisim.persistence +from romanisim import log +import roman_datamodels.maker_utils as maker_utils + +# Define centermost SCA for PSFs +CENTER_SCA = 2 def add_objects_to_l3(l3_mos, source_cat, rng=None, seed=None): @@ -35,46 +53,398 @@ def add_objects_to_l3(l3_mos, source_cat, rng=None, seed=None): filter_name = l3_mos.meta.basic.optical_element # Generate WCS - twcs = wcs.get_mosaic_wcs(l3_mos.meta) + twcs = romanisim.wcs.get_mosaic_wcs(l3_mos.meta, shape=l3_mos.data.shape) # Create PSF - l3_psf = psf.make_psf(filter_name=filter_name, sca=2, chromatic=False, webbpsf=True) + l3_psf = romanisim.psf.make_psf(filter_name=filter_name, sca=CENTER_SCA, chromatic=False, webbpsf=True) # Generate x,y positions for sources coords = np.array([[o.sky_pos.ra.rad, o.sky_pos.dec.rad] for o in source_cat]) sourcecountsall = galsim.ImageF(l3_mos.data.shape[0], l3_mos.data.shape[1], wcs=twcs, xmin=0, ymin=0) - xpos, ypos = sourcecountsall.wcs._xy(coords[:, 0], coords[:, 1]) + # xpos, ypos = sourcecountsall.wcs._xy(coords[:, 0], coords[:, 1]) + xpos, ypos = sourcecountsall.wcs.radecToxy(coords[:, 0], coords[:, 1], 'rad') xpos_idx = [round(x) for x in xpos] ypos_idx = [round(y) for y in ypos] + if ((min(xpos_idx) < 0) or (min(ypos_idx)) or (max(xpos_idx) > l3_mos.data.shape[0]) + or (max(ypos_idx) > l3_mos.data.shape[1])): + log.error(f"A source is out of bounds! " + f"Source min x,y = {min(xpos_idx)},{min(ypos_idx)}, max = {max(xpos_idx)},{max(ypos_idx)} " + f"Image min x,y = {0},{0}, max = {l3_mos.data.shape[0]},{l3_mos.data.shape[1]}") + # Create overall scaling factor map - Ct_all = (l3_mos.data.value / l3_mos.var_poisson) + # Ct_all = (l3_mos.data.value / l3_mos.var_poisson) + Ct_all = np.divide(l3_mos.data.value, l3_mos.var_poisson.value, + out=np.ones(l3_mos.data.shape), where=l3_mos.var_poisson.value != 0) # Cycle over sources and add them to the mosaic for idx, (x, y) in enumerate(zip(xpos_idx, ypos_idx)): # Set scaling factor for injected sources # Flux / sigma_p^2 - Ct = (l3_mos.data[x][y].value / l3_mos.var_poisson[x][y].value) + if l3_mos.var_poisson[x][y].value != 0: + Ct = math.fabs(l3_mos.data[x][y].value / l3_mos.var_poisson[x][y].value) + # elif l3_mos.data[x][y].value != 0: + # Ct = math.fabs(l3_mos.data[x][y].value) + else: + continue # Create empty postage stamp galsim source image sourcecounts = galsim.ImageF(l3_mos.data.shape[0], l3_mos.data.shape[1], wcs=twcs, xmin=0, ymin=0) # Simulate source postage stamp - image.add_objects_to_image(sourcecounts, [source_cat[idx]], xpos=[xpos[idx]], ypos=[ypos[idx]], - psf=l3_psf, flux_to_counts_factor=Ct, bandpass=[filter_name], - filter_name=filter_name, rng=rng) + romanisim.image.add_objects_to_image(sourcecounts, [source_cat[idx]], xpos=[xpos[idx]], + ypos=[ypos[idx]], psf=l3_psf, flux_to_counts_factor=Ct, + bandpass=[filter_name], filter_name=filter_name, rng=rng) # Scale the source image back by its flux ratios sourcecounts /= Ct # Add sources to the original mosaic data array - l3_mos.data = (l3_mos.data.value + sourcecounts.array) * l3_mos.data.unit + # l3_mos.data = (l3_mos.data.value + sourcecounts.array) * l3_mos.data.unit + l3_mos.data = (l3_mos.data.value + np.swapaxes(sourcecounts.array, 0, 1)) * l3_mos.data.unit # Note for the future - other noise sources (read and flat) need to be implemented # Set new poisson variance - l3_mos.var_poisson = l3_mos.data.value / Ct_all + # l3_mos.var_poisson = l3_mos.data.value / Ct_all + l3_mos.var_poisson = np.divide(l3_mos.data.value, Ct_all, out=l3_mos.var_poisson.value, + where=(Ct_all != 0)) * l3_mos.var_poisson.unit # l3_mos is updated in place, so no return return None + + +def make_l3(image, context, metadata, exptimes, rng=None, seed=None, + saturation=None, unit_factor=1): + """TBD + """ + + log.info('Apportioning counts to mosaic...') + # resultants, dq = apportion_counts_to_resultants( + # counts.array, tij, inv_linearity=inv_linearity, crparam=crparam, + # persistence=persistence, tstart=tstart, + # rng=rng, seed=seed) + + # roman.addReciprocityFailure(resultants_object) + + # code from apportion_counts_to_resultants + + mosaic = np.swapaxes(image.array, 0, 1) + + # Set rng for creating readnoise, flat noise + if rng is None and seed is None: + seed = 46 + log.warning( + 'No RNG set, constructing a new default RNG from default seed.') + if rng is None: + rng = galsim.GaussianDeviate(seed) + + # Need to add read noise, poisson noise, and flat noise + + # Simulate read and flat as random noise + # if read_noise is not None and not isinstance(read_noise, u.Quantity): + # read_noise = read_noise * u.DN + # log.warning('Making up units for read noise.') + # resultants are now in counts. + # read noise is in counts. + # log.info('Adding read noise...') + # resultants = add_read_noise_to_resultants( + # resultants, tij, rng=rng, seed=seed, + # read_noise=read_noise) + + # Improve this for proper 3D contexts + # ctx_2d = np.reshape(context[-2:], context.shape[-2:]) + # mosaic = np.divide(counts_arr, context[-2:], where=context[-2:]!=0, out=np.zeros(counts_arr.shape, dtype='f8')) + # mosaic = np.divide(counts_arr, ctx_2d, where=ctx_2d!=0, out=np.zeros(counts_arr.shape, dtype='f8')) + + # Mosaics may suffer saturation? If so rework the following. + # if saturation is None: + # saturation = parameters.reference_data['saturation'] + # mosaic = np.clip(mosaic, 0, 2 * 10**9).astype('i4') + # this maybe should be better applied at read time? + # it's not actually clear to me what the right thing to do + # is in detail. + # mosaic = np.clip(mosaic, 0 * u.DN, saturation) + # mosaic = np.clip(mosaic, 0, saturation.value) + # m = mosaic >= saturation.value + # dq = np.zeros(counts_arr.shape, dtype=np.uint32) + # dq[m] |= parameters.dqbits['saturated'] + + # Set mosaic to be a mosaic model + mosaic_mdl = maker_utils.mk_level3_mosaic(shape=mosaic.shape, meta=metadata) + + # Extract exposure time for each pixel + exptimes_pix = util.decode_context_times(context, exptimes) + + # Set data + mosaic_mdl.data = mosaic * unit_factor + mosaic_mdl.data = np.divide(mosaic_mdl.data.value, exptimes_pix, + out=np.zeros(mosaic_mdl.data.shape), + where=exptimes_pix != 0) * mosaic_mdl.data.unit + + # Context + # Binary index of images that contribute to the pixel + # Defined by geometry.. if not, all non zero = 1. + mosaic_mdl.context = context + + # Poisson noise + mosaic_mdl.var_poisson = unit_factor**2 / exptimes_pix**2 + + # Weight + # Use exptime weight + # --- + # DQ used if set above for saturation? + # bitvalue = interpret_bit_flags( + # bitvalue, flag_name_map={dq.name: dq.value for dq in pixel} ) + # if bitvalue is None: + # return np.ones(dqarr.shape, dtype=np.uint8) + # return np.logical_not(np.bitwise_and(dqarr, ~bitvalue)).astype(np.uint8) + # --- + # dqmask = build_mask(model.dq, good_bits) + # elif weight_type == "exptime": + # exptime = model.meta.exposure.exposure_time + # result = exptime * dqmask + + # Fill everything else: + # err, weight, var_rnoise, var_flat + + return mosaic_mdl + + +def generate_mosaic_geometry(): + """Generate a geometry map (context) for a mosaic dither pattern / set of pointings and roll angles + TBD + """ + pass + + +def simulate(metadata, cat, exptimes, context=None, + usecrds=True, webbpsf=True, seed=None, rng=None, + psf_keywords=dict(), **kwargs + ): + """TBD + """ + + 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.') + + # Create metadata object + meta = maker_utils.mk_mosaic_meta() + meta['wcs'] = None + + for key in parameters.default_mosaic_parameters_dictionary.keys(): + meta[key].update(parameters.default_mosaic_parameters_dictionary[key]) + + for key in metadata.keys(): + if key not in meta: + meta[key] = metadata[key] + # elif isinstance(meta[key], dict): + else: + meta[key].update(metadata[key]) + + # May need an equivalent of this this for L3 + # util.add_more_metadata(meta) + + filter_name = metadata['basic']['optical_element'] + + if rng is None and seed is None: + seed = 43 + log.warning( + 'No RNG set, constructing a new default RNG from default seed.') + if rng is None: + rng = galsim.UniformDeviate(seed) + + if context is None: + # Create geometry from the object list + twcs = romanisim.wcs.get_mosaic_wcs(meta) + + coords = np.array([[o.sky_pos.ra.rad, o.sky_pos.dec.rad] + for o in cat]) + + allx, ally = twcs.radecToxy(coords[:, 0], coords[:, 1], 'rad') + + # Obtain the sample extremums + xmin = min(allx) + xmax = max(allx) + ymin = min(ally) + ymax = max(ally) + + # Obtain WCS center + xcen, ycen = twcs.radecToxy(twcs.center.ra, twcs.center.dec, 'rad') + + # Determine maximum extremums from WCS center + xdiff = max([math.ceil(xmax - xcen), math.ceil(xcen - xmin)]) + ydiff = max([math.ceil(ymax - ycen), math.ceil(ycen - ymin)]) + + # Create context map preserving WCS center + context = np.ones((1, 2 * xdiff, 2 * ydiff), dtype=np.uint32) + + log.info('Simulating filter {0}...'.format(filter_name)) + mosaic, simcatobj = simulate_cps( + meta, cat, context, exptimes, rng=rng, usecrds=usecrds, + webbpsf=webbpsf, psf_keywords=psf_keywords) + + extras = {} + + # if reffiles: + # extras["simulate_reffiles"] = {} + # for key, value in reffiles.items(): + # extras["simulate_reffiles"][key] = value + + extras['simcatobj'] = simcatobj + + # TODO: Double check that this is a gwcs object + # extras['wcs'] = wcs.convert_wcs_to_gwcs(mosaic.wcs) + # extras['wcs'] = mosaic.wcs + + log.info('Simulation complete.') + # return mosaic, extras + return mosaic, extras + + +def simulate_cps(metadata, objlist, context, exptimes, + rng=None, seed=None, webbpsf=True, usecrds=False, + psf_keywords=dict()): + """TBD + """ + + if rng is None and seed is None: + seed = 43 + log.warning( + 'No RNG set, constructing a new default RNG from default seed.') + if rng is None: + rng = galsim.UniformDeviate(seed) + + filter_name = metadata['basic']['optical_element'] + + date = metadata['basic']['time_mean_mjd'] + + if not isinstance(date, astropy.time.Time): + date = astropy.time.Time(date, format='mjd') + + galsim_filter_name = romanisim.bandpass.roman2galsim_bandpass[filter_name] + bandpass = roman.getBandpasses(AB_zeropoint=True)[galsim_filter_name] + + # Generate WCS + moswcs = romanisim.wcs.get_mosaic_wcs(metadata, shape=context.shape[-2:]) + + # Is this needed? + chromatic = False + if (len(objlist) > 0 + and not isinstance(objlist, table.Table) # this case is always gray + and objlist[0].profile.spectral): + chromatic = True + + # With context, is possible that individual PSFs could be properly assigned? + # Would be a pain to implement + psf = romanisim.psf.make_psf(filter_name=filter_name, wcs=moswcs, sca=CENTER_SCA, + chromatic=chromatic, webbpsf=webbpsf, + variable=True, **psf_keywords) + + # Create initial galsim image + image = galsim.ImageF(context.shape[-2], context.shape[-1], wcs=moswcs, xmin=0, ymin=0) + + # Examine this in more detail to ensure it is correct + # Create base sky level + mos_cent_pos = moswcs.toWorld(image.true_center) + sky_level = roman.getSkyLevel(bandpass, world_pos=mos_cent_pos, + date=date.datetime, exptime=1) + sky_level *= (1.0 + roman.stray_light_fraction) + # sky_mosaic = image * 0 + # moswcs.makeSkyImage(sky_mosaic, sky_level) + # sky_mosaic += roman.thermal_backgrounds[galsim_filter_name] + abflux = romanisim.bandpass.get_abflux(filter_name) + + # Obtain physical unit conversion factor + unit_factor = parameters.reference_data['photom'][filter_name] + + # Maybe return a mosaic object as well? + mosaic, simcatobj = simulate_cps_generic( + image, metadata, context, exptimes, objlist=objlist, psf=psf, zpflux=abflux, sky=sky_level, + wcs=moswcs, rng=rng, seed=seed, unit_factor=unit_factor) + + return mosaic, simcatobj + + +def simulate_cps_generic(image, metadata, context, exptimes, objlist=None, psf=None, + zpflux=None, wcs=None, xpos=None, ypos=None, sky=None, + flat=None, rng=None, seed=None, unit_factor=1): + """TBD + """ + if rng is None and seed is None: + seed = 144 + log.warning( + 'No RNG set, constructing a new default RNG from default seed.') + if rng is None: + rng = galsim.UniformDeviate(seed) + + # Objlist wcs check + if (objlist is not None and len(objlist) > 0 + and image.wcs is None and (xpos is None or ypos is None)): + raise ValueError('xpos and ypos must be set if rendering objects ' + 'without a WCS.') + if objlist is None: + objlist = [] + if len(objlist) > 0 and xpos is None: + if isinstance(objlist, table.Table): + xpos, ypos = image.wcs.radecToxy( + np.radians(objlist['ra']), np.radians(objlist['dec']), 'rad') + else: + coord = np.array([[o.sky_pos.ra.rad, o.sky_pos.dec.rad] + for o in objlist]) + xpos, ypos = image.wcs.radecToxy(coord[:, 0], coord[:, 1], 'rad') + # use private vectorized transformation + if xpos is not None: + xpos = np.array(xpos) + if ypos is not None: + ypos = np.array(ypos) + + if len(objlist) > 0 and psf is None: + raise ValueError('Must provide a PSF if you want to render objects.') + + if flat is None: + flat = 1.0 + + # # for some reason, galsim doesn't like multiplying an SED by 1, but it's + # # okay with multiplying an SED by 1.0, so we cast to float. + 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. + + chromatic = False + if len(objlist) > 0 and objlist[0].profile.spectral: + chromatic = True + + # add Poisson noise if we made a noiseless, not-photon-shooting + # image. + poisson_noise = galsim.PoissonNoise(rng, sky_level=sky) + if not chromatic: + image.addNoise(poisson_noise) + + # Add Objects requires a mosaic image + # Maybe violates the notion of a "generic" count simulator? + mosaic_mdl = make_l3(image, context, metadata, exptimes, unit_factor=unit_factor) + + if wcs is not None: + mosaic_mdl['wcs'] = wcs + + # Add objects to mosaic model + objinfo = add_objects_to_l3( + mosaic_mdl, objlist, rng=rng) + objinfo = np.zeros( + len(objlist), + dtype=[('x', 'f4'), ('y', 'f4'), ('counts', 'f4'), ('time', 'f4')]) + if len(objlist) > 0: + objinfo['x'] = xpos + objinfo['y'] = ypos + + return mosaic_mdl, objinfo diff --git a/romanisim/parameters.py b/romanisim/parameters.py index 221a14f2..67ba1df9 100644 --- a/romanisim/parameters.py +++ b/romanisim/parameters.py @@ -84,6 +84,15 @@ "linearity": None, "readnoise": 5.0 * u.DN, "saturation": 55000 * u.DN, + "photom": {"F062": 0.28798833 * u.MJy / u.sr, + "F087": 0.3911463 * u.MJy / u.sr, + "F106": 0.37039083 * u.MJy / u.sr, + "F129": 0.36495176 * u.MJy / u.sr, + "F146": 0.11750617 * u.MJy / u.sr, + "F158": 0.35352992 * u.MJy / u.sr, + "F184": 0.53612772 * u.MJy / u.sr, + "F213": 0.55763922 * u.MJy / u.sr, + } } nborder = 4 # number of border pixels used for reference pixels. diff --git a/romanisim/tests/test_l3.py b/romanisim/tests/test_l3.py new file mode 100644 index 00000000..15b03ae2 --- /dev/null +++ b/romanisim/tests/test_l3.py @@ -0,0 +1,160 @@ +"""Unit tests for mosaic module. + +""" + +import os +import copy +import numpy as np +import galsim +from galsim import roman +from romanisim import image, parameters, catalog, psf, util, wcs, persistence, ramp, l1, l3 +from astropy.coordinates import SkyCoord +from astropy import units as u +from astropy.time import Time +from astropy import table +import asdf +import webbpsf +from astropy.modeling.functional_models import Sersic2D +import pytest +from metrics_logger.decorators import metrics_logger +from romanisim import log +from roman_datamodels.stnode import WfiScienceRaw, WfiImage +import roman_datamodels.maker_utils as maker_utils +import romanisim.bandpass + + +@metrics_logger("DMS232") +@pytest.mark.soctests +def test_inject_sources_into_mosaic1(): + """Inject sources into a mosaic. + """ + + # Set constants and metadata + galsim.roman.n_pix = 200 + rng_seed = 42 + metadata = copy.deepcopy(parameters.default_mosaic_parameters_dictionary) + metadata['basic']['optical_element'] = 'F158' + + # Create WCS + twcs = wcs.get_mosaic_wcs(metadata, shape=(galsim.roman.n_pix, galsim.roman.n_pix)) + + # Create initial Level 3 mosaic + + # Create Four-quadrant pattern of gaussian noise, centered around one + # Each quadrant's gaussian noise scales like total exposure time + # (total files contributed to each quadrant) + + # Create gaussian noise generators + g1 = galsim.GaussianDeviate(rng_seed, mean=1.0, sigma=0.01) + g2 = galsim.GaussianDeviate(rng_seed, mean=1.0, sigma=0.02) + g3 = galsim.GaussianDeviate(rng_seed, mean=1.0, sigma=0.05) + g4 = galsim.GaussianDeviate(rng_seed, mean=1.0, sigma=0.1) + + # Create level 3 mosaic model + l3_mos = maker_utils.mk_level3_mosaic(shape=(galsim.roman.n_pix, galsim.roman.n_pix)) + + # Update metadata in the l3 model + for key in metadata.keys(): + if key in l3_mos.meta: + l3_mos.meta[key].update(metadata[key]) + + # Populate the mosaic data array with gaussian noise from generators + g1.generate(l3_mos.data.value[0:100, 0:100]) + g2.generate(l3_mos.data.value[0:100, 100:200]) + g3.generate(l3_mos.data.value[100:200, 0:100]) + g4.generate(l3_mos.data.value[100:200, 100:200]) + + # Define Poisson Noise of mosaic + l3_mos.var_poisson.value[0:100, 0:100] = 0.01**2 + l3_mos.var_poisson.value[0:100, 100:200] = 0.02**2 + l3_mos.var_poisson.value[100:200, 0:100] = 0.05**2 + l3_mos.var_poisson.value[100:200, 100:200] = 0.1**2 + + # Create normalized psf source catalog (same source in each quadrant) + sc_dict = {"ra": 4 * [0.0], "dec": 4 * [0.0], "type": 4 * ["PSF"], "n": 4 * [-1.0], + "half_light_radius": 4 * [0.0], "pa": 4 * [0.0], "ba": 4 * [1.0], "F158": 4 * [1.0]} + sc_table = table.Table(sc_dict) + + xpos, ypos = 50, 50 + sc_table["ra"][0], sc_table["dec"][0] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value + xpos, ypos = 50, 150 + sc_table['ra'][1], sc_table['dec'][1] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value + xpos, ypos = 150, 50 + sc_table['ra'][2], sc_table['dec'][2] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value + xpos, ypos = 150, 150 + sc_table['ra'][3], sc_table['dec'][3] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value + + source_cat = catalog.table_to_catalog(sc_table, ["F158"]) + + # Copy original Mosaic before adding sources + l3_mos_orig = l3_mos.copy() + + # Add source_cat objects to mosaic + l3.add_objects_to_l3(l3_mos, source_cat, seed=rng_seed) + + # Ensure that every data pixel value has increased or + # remained the same with the new sources injected + assert np.all(l3_mos.data.value >= l3_mos_orig.data.value) + + # Ensure that every pixel's poisson variance has increased or + # remained the same with the new sources injected + # Numpy isclose is needed to determine equality, due to float precision issues + close_mask = np.isclose(l3_mos.var_poisson.value, l3_mos_orig.var_poisson.value, rtol=1e-06) + assert np.all(l3_mos.var_poisson.value[~close_mask] > l3_mos_orig.var_poisson.value[~close_mask]) + + # Create log entry and artifacts + log.info('DMS232 successfully injected sources into a mosaic at points (50,50), (50,150), (150,50), (150,150).') + + artifactdir = os.environ.get('TEST_ARTIFACT_DIR', None) + if artifactdir is not None: + af = asdf.AsdfFile() + af.tree = {'l3_mos': l3_mos, + 'l3_mos_orig': l3_mos_orig, + 'source_cat_table': sc_table, + } + af.write_to(os.path.join(artifactdir, 'dms232.asdf')) + + +def test_sim_mosaic1(): + """Simulating mosaic from catalog file. + """ + + ra_ref = 1.0 + dec_ref = 4.0 + + metadata = copy.deepcopy(parameters.default_mosaic_parameters_dictionary) + metadata['basic']['optical_element'] = 'F158' + metadata['wcsinfo']['ra_ref'] = ra_ref + metadata['wcsinfo']['dec_ref'] = dec_ref + metadata['wcsinfo']['pixel_scale'] = 0.11 + metadata['wcsinfo']['pixel_scale_local'] = 0.11 + metadata['wcsinfo']['v2_ref'] = 0 + metadata['wcsinfo']['v3_ref'] = 0 + + exptimes = [600] + + cen = SkyCoord(ra=ra_ref * u.deg, dec=dec_ref * u.deg) + cat = catalog.make_dummy_table_catalog(cen, radius=0.01, nobj=100) + cat['F158'] = cat['F158'] * 10e10 + + source_cat = catalog.table_to_catalog(cat, ["F158"]) + + mosaic, extras = l3.simulate(metadata, source_cat, exptimes) + + import plotly.express as px + + fig1 = px.imshow(mosaic.data.value, title='Mosaic Data', labels={'color': 'MJy / sr'}) + fig1.show() + + fig2 = px.imshow(mosaic.context[-2:].reshape(mosaic.context.shape[-2:]), title='Mosaic Context', labels={'color': 'File Number'}) + fig2.show() + + pos_vals = mosaic.data.value.copy() + pos_vals[pos_vals <= 0] = 0.00000000001 + + fig3 = px.imshow(np.log(pos_vals), title='Mosaic Data (log)', labels={'color': 'MJy / sr'}) + fig3.show() + +# TBD: Test with more complex context + +# TBD: Test of geometry construction diff --git a/romanisim/util.py b/romanisim/util.py index 82a8f0f1..9bf57dc2 100644 --- a/romanisim/util.py +++ b/romanisim/util.py @@ -298,3 +298,87 @@ def sample_king_distances(rc, rt, npts, rng=None): cdf /= cdf[-1] radii = np.interp(rr, cdf, x) return radii + + +def decode_context_times(context, exptimes): + """ + Get 0-based indices of input images that contributed to (resampled) + output pixel with coordinates ``x`` and ``y``. + + Parameters + ---------- + context: numpy.ndarray + A 3D `~numpy.ndarray` of integral data type. + + exptimes: list of floats + Exposure times for each component image. + + + Returns + ------- + total_exptimes: numpy.ndarray + A 2D array of total exposure time for each pixel. + """ + + if context.ndim != 3: + raise ValueError("'context' must be a 3D array.") + + """ + Context decoding example: + An example context array for an output image of array shape ``(5, 6)`` + obtained by resampling 80 input images. + + .. code-block:: none + + con = np.array( + [[[0, 0, 0, 0, 0, 0], + [0, 0, 0, 36196864, 0, 0], + [0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0], + [0, 0, 537920000, 0, 0, 0]], + [[0, 0, 0, 0, 0, 0,], + [0, 0, 0, 67125536, 0, 0], + [0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0], + [0, 0, 163856, 0, 0, 0]], + [[0, 0, 0, 0, 0, 0], + [0, 0, 0, 8203, 0, 0], + [0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0], + [0, 0, 32865, 0, 0, 0]]], + dtype=np.int32 + ) + decode_context(con, [3, 2], [1, 4]) + [array([ 9, 12, 14, 19, 21, 25, 37, 40, 46, 58, 64, 65, 67, 77]), + array([ 9, 20, 29, 36, 47, 49, 64, 69, 70, 79])] + """ + + nbits = 8 * context.dtype.itemsize + + total_exptimes = np.zeros(context.shape[1:]) + + for x in range(total_exptimes.shape[0]): + for y in range(total_exptimes.shape[1]): + files = [v & (1 << k) for v in context[:, x, y] for k in range(nbits)] + tot_time = 0 + files = [file for file in files if (file != 0)] + + for im_idx in files: + tot_time += exptimes[im_idx - 1] + + total_exptimes[x][y] = tot_time + + def sum_times(x): + tot_time = 0 + files = [x & (1 << k) for k in range(nbits)] + files = [file for file in files if (file != 0)] + for im_idx in files: + tot_time += exptimes[im_idx - 1] + return tot_time + + vectorized_sum_times = np.vectorize(sum_times) + + total_exptimes = vectorized_sum_times(context[:,]) + total_exptimes = total_exptimes.reshape(total_exptimes.shape[1:]) + + return total_exptimes diff --git a/romanisim/wcs.py b/romanisim/wcs.py index 04830824..62a09870 100644 --- a/romanisim/wcs.py +++ b/romanisim/wcs.py @@ -17,6 +17,7 @@ """ import warnings +import math import numpy as np import astropy.coordinates from astropy import units as u @@ -408,7 +409,7 @@ def convert_wcs_to_gwcs(wcs): return wcs_from_fits_header(wcs.header.header) -def get_mosaic_wcs(mosaic): +def get_mosaic_wcs(mosaic, shape=None): """Get a WCS object for a given sca or set of CRDS parameters. Parameters @@ -423,7 +424,10 @@ def get_mosaic_wcs(mosaic): # If sent a dictionary, create a temporary model for data interface if (type(mosaic) is not roman_datamodels.datamodels.MosaicModel): - mosaic_node = maker_utils.mk_level3_mosaic() + if shape is None: + mosaic_node = maker_utils.mk_level3_mosaic() + else: + mosaic_node = maker_utils.mk_level3_mosaic(shape=shape) for key in mosaic.keys(): if isinstance(mosaic[key], dict): mosaic_node['meta'][key].update(mosaic[key]) @@ -437,13 +441,17 @@ def get_mosaic_wcs(mosaic): mosaic_node.meta.wcsinfo.ra_ref * u.deg, mosaic_node.meta.wcsinfo.dec_ref * u.deg) + if shape is None: + shape = (mosaic_node.data.shape[0], + mosaic_node.data.shape[1]) + # Create a tangent plane WCS for the mosaic # The affine parameters below should be reviewed and updated affine = galsim.AffineTransform( - 0.1, 0, 0, 0.1, origin=galsim.PositionI(mosaic_node.data.shape[0] / 2.0, - mosaic_node.data.shape[1] / 2.0,), + 0.1, 0, 0, 0.1, origin=galsim.PositionI(math.ceil(shape[0] / 2.0), math.ceil(shape[1] / 2.0)), world_origin=galsim.PositionD(0, 0)) wcs = galsim.TanWCS(affine, util.celestialcoord(world_pos)) + # wcs = GWCS(wcs) return wcs From 9e3e98b510431befd52f6e701cb5e6bfee486cee Mon Sep 17 00:00:00 2001 From: PaulHuwe Date: Mon, 24 Jun 2024 13:36:59 -0400 Subject: [PATCH 2/9] Moved L3 tests to test_l3 --- romanisim/tests/test_image.py | 92 ----------------------------------- romanisim/tests/test_l3.py | 4 +- 2 files changed, 2 insertions(+), 94 deletions(-) diff --git a/romanisim/tests/test_image.py b/romanisim/tests/test_image.py index 554b2227..c56378f2 100644 --- a/romanisim/tests/test_image.py +++ b/romanisim/tests/test_image.py @@ -755,95 +755,3 @@ def test_inject_source_into_image(): 'flux': flux, 'tij': tij} af.write_to(os.path.join(artifactdir, 'dms231.asdf')) - - -@metrics_logger("DMS232") -@pytest.mark.soctests -def test_inject_sources_into_mosaic(): - """Inject sources into a mosaic. - """ - - # Set constants and metadata - galsim.roman.n_pix = 200 - rng_seed = 42 - metadata = copy.deepcopy(parameters.default_mosaic_parameters_dictionary) - metadata['basic']['optical_element'] = 'F158' - - # Create WCS - twcs = wcs.get_mosaic_wcs(metadata) - - # Create initial Level 3 mosaic - - # Create Four-quadrant pattern of gaussian noise, centered around one - # Each quadrant's gaussian noise scales like total exposure time - # (total files contributed to each quadrant) - - # Create gaussian noise generators - g1 = galsim.GaussianDeviate(rng_seed, mean=1.0, sigma=0.01) - g2 = galsim.GaussianDeviate(rng_seed, mean=1.0, sigma=0.02) - g3 = galsim.GaussianDeviate(rng_seed, mean=1.0, sigma=0.05) - g4 = galsim.GaussianDeviate(rng_seed, mean=1.0, sigma=0.1) - - # Create level 3 mosaic model - l3_mos = maker_utils.mk_level3_mosaic(shape=(galsim.roman.n_pix, galsim.roman.n_pix)) - - # Update metadata in the l3 model - for key in metadata.keys(): - if key in l3_mos.meta: - l3_mos.meta[key].update(metadata[key]) - - # Populate the mosaic data array with gaussian noise from generators - g1.generate(l3_mos.data.value[0:100, 0:100]) - g2.generate(l3_mos.data.value[0:100, 100:200]) - g3.generate(l3_mos.data.value[100:200, 0:100]) - g4.generate(l3_mos.data.value[100:200, 100:200]) - - # Define Poisson Noise of mosaic - l3_mos.var_poisson.value[0:100, 0:100] = 0.01**2 - l3_mos.var_poisson.value[0:100, 100:200] = 0.02**2 - l3_mos.var_poisson.value[100:200, 0:100] = 0.05**2 - l3_mos.var_poisson.value[100:200, 100:200] = 0.1**2 - - # Create normalized psf source catalog (same source in each quadrant) - sc_dict = {"ra": 4 * [0.0], "dec": 4 * [0.0], "type": 4 * ["PSF"], "n": 4 * [-1.0], - "half_light_radius": 4 * [0.0], "pa": 4 * [0.0], "ba": 4 * [1.0], "F158": 4 * [1.0]} - sc_table = table.Table(sc_dict) - - xpos, ypos = 50, 50 - sc_table["ra"][0], sc_table["dec"][0] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value - xpos, ypos = 50, 150 - sc_table['ra'][1], sc_table['dec'][1] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value - xpos, ypos = 150, 50 - sc_table['ra'][2], sc_table['dec'][2] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value - xpos, ypos = 150, 150 - sc_table['ra'][3], sc_table['dec'][3] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value - - source_cat = catalog.table_to_catalog(sc_table, ["F158"]) - - # Copy original Mosaic before adding sources - l3_mos_orig = l3_mos.copy() - - # Add source_cat objects to mosaic - l3.add_objects_to_l3(l3_mos, source_cat, seed=rng_seed) - - # Ensure that every data pixel value has increased or - # remained the same with the new sources injected - assert np.all(l3_mos.data.value >= l3_mos_orig.data.value) - - # Ensure that every pixel's poisson variance has increased or - # remained the same with the new sources injected - # Numpy isclose is needed to determine equality, due to float precision issues - close_mask = np.isclose(l3_mos.var_poisson.value, l3_mos_orig.var_poisson.value, rtol=1e-06) - assert np.all(l3_mos.var_poisson.value[~close_mask] > l3_mos_orig.var_poisson.value[~close_mask]) - - # Create log entry and artifacts - log.info('DMS232 successfully injected sources into a mosiac at points (50,50), (50,150), (150,50), (150,150).') - - artifactdir = os.environ.get('TEST_ARTIFACT_DIR', None) - if artifactdir is not None: - af = asdf.AsdfFile() - af.tree = {'l3_mos': l3_mos, - 'l3_mos_orig': l3_mos_orig, - 'source_cat_table': sc_table, - } - af.write_to(os.path.join(artifactdir, 'dms232.asdf')) diff --git a/romanisim/tests/test_l3.py b/romanisim/tests/test_l3.py index 15b03ae2..6cc580ef 100644 --- a/romanisim/tests/test_l3.py +++ b/romanisim/tests/test_l3.py @@ -25,7 +25,7 @@ @metrics_logger("DMS232") @pytest.mark.soctests -def test_inject_sources_into_mosaic1(): +def test_inject_sources_into_mosaic(): """Inject sources into a mosaic. """ @@ -115,7 +115,7 @@ def test_inject_sources_into_mosaic1(): af.write_to(os.path.join(artifactdir, 'dms232.asdf')) -def test_sim_mosaic1(): +def test_sim_mosaic(): """Simulating mosaic from catalog file. """ From 0d3e83c5579820d5fd316c0231666f9db0e2eecd Mon Sep 17 00:00:00 2001 From: PaulHuwe Date: Sat, 6 Jul 2024 14:32:42 -0400 Subject: [PATCH 3/9] Save --- romanisim/l3.py | 44 +++++++++++++++++++++----------------------- 1 file changed, 21 insertions(+), 23 deletions(-) diff --git a/romanisim/l3.py b/romanisim/l3.py index db5a5af0..84faaed2 100644 --- a/romanisim/l3.py +++ b/romanisim/l3.py @@ -26,7 +26,7 @@ CENTER_SCA = 2 -def add_objects_to_l3(l3_mos, source_cat, rng=None, seed=None): +def add_objects_to_l3(l3_mos, xpos=None, ypos=None, coords=None, coords_unit='rad', wcs=None, psf=None, rng=None, seed=None): """Add objects to a Level 3 mosaic Parameters @@ -42,28 +42,25 @@ def add_objects_to_l3(l3_mos, source_cat, rng=None, seed=None): l3_mos is updated in place """ - if rng is None and seed is None: - seed = 143 - log.warning( - 'No RNG set, constructing a new default RNG from default seed.') - if rng is None: - rng = galsim.UniformDeviate(seed) - # Obtain optical element filter_name = l3_mos.meta.basic.optical_element - # Generate WCS - twcs = romanisim.wcs.get_mosaic_wcs(l3_mos.meta, shape=l3_mos.data.shape) + if wcs is None: + # Generate WCS + wcs = romanisim.wcs.get_mosaic_wcs(l3_mos.meta, shape=l3_mos.data.shape) + + if psf is None: + # Create PSF + psf = romanisim.psf.make_psf(filter_name=filter_name, sca=CENTER_SCA, chromatic=False, webbpsf=True) - # Create PSF - l3_psf = romanisim.psf.make_psf(filter_name=filter_name, sca=CENTER_SCA, chromatic=False, webbpsf=True) + if coords is None: + # Generate x,y positions for sources + sourcecountsall = galsim.ImageF(l3_mos.data.shape[0], l3_mos.data.shape[1], wcs=wcs, xmin=0, ymin=0) + xpos, ypos = sourcecountsall.wcs.radecToxy(coords[:, 0], coords[:, 1], coords_unit) + elif all(pos is not None for pos in [xpos,ypos]): + log.error('No (xpos and ypos) or coords set. Adding zero objects.') + return - # Generate x,y positions for sources - coords = np.array([[o.sky_pos.ra.rad, o.sky_pos.dec.rad] - for o in source_cat]) - sourcecountsall = galsim.ImageF(l3_mos.data.shape[0], l3_mos.data.shape[1], wcs=twcs, xmin=0, ymin=0) - # xpos, ypos = sourcecountsall.wcs._xy(coords[:, 0], coords[:, 1]) - xpos, ypos = sourcecountsall.wcs.radecToxy(coords[:, 0], coords[:, 1], 'rad') xpos_idx = [round(x) for x in xpos] ypos_idx = [round(y) for y in ypos] @@ -90,19 +87,19 @@ def add_objects_to_l3(l3_mos, source_cat, rng=None, seed=None): continue # Create empty postage stamp galsim source image - sourcecounts = galsim.ImageF(l3_mos.data.shape[0], l3_mos.data.shape[1], wcs=twcs, xmin=0, ymin=0) + sourcecounts = galsim.ImageF(l3_mos.data.shape[0], l3_mos.data.shape[1], wcs=wcs, xmin=0, ymin=0) # Simulate source postage stamp romanisim.image.add_objects_to_image(sourcecounts, [source_cat[idx]], xpos=[xpos[idx]], - ypos=[ypos[idx]], psf=l3_psf, flux_to_counts_factor=Ct, + ypos=[ypos[idx]], psf=psf, flux_to_counts_factor=Ct, bandpass=[filter_name], filter_name=filter_name, rng=rng) # Scale the source image back by its flux ratios sourcecounts /= Ct # Add sources to the original mosaic data array - # l3_mos.data = (l3_mos.data.value + sourcecounts.array) * l3_mos.data.unit - l3_mos.data = (l3_mos.data.value + np.swapaxes(sourcecounts.array, 0, 1)) * l3_mos.data.unit + l3_mos.data = (l3_mos.data.value + sourcecounts.array) * l3_mos.data.unit + # l3_mos.data = (l3_mos.data.value + np.swapaxes(sourcecounts.array, 0, 1)) * l3_mos.data.unit # Note for the future - other noise sources (read and flat) need to be implemented @@ -130,7 +127,8 @@ def make_l3(image, context, metadata, exptimes, rng=None, seed=None, # code from apportion_counts_to_resultants - mosaic = np.swapaxes(image.array, 0, 1) + # mosaic = np.swapaxes(image.array, 0, 1) + mosaic = image.array.copy() # Set rng for creating readnoise, flat noise if rng is None and seed is None: From e1de7e34a0a1749b2cd97c0fcb3df99c32adecd5 Mon Sep 17 00:00:00 2001 From: PaulHuwe Date: Sun, 7 Jul 2024 15:43:40 -0400 Subject: [PATCH 4/9] Save --- romanisim/image.py | 132 +++++++++++- romanisim/l3.py | 413 +++++-------------------------------- romanisim/tests/test_l3.py | 67 ++---- 3 files changed, 196 insertions(+), 416 deletions(-) diff --git a/romanisim/image.py b/romanisim/image.py index 0727f3e5..d70c03cd 100644 --- a/romanisim/image.py +++ b/romanisim/image.py @@ -210,7 +210,8 @@ def trim_objlist(objlist, image): def add_objects_to_image(image, objlist, xpos, ypos, psf, - flux_to_counts_factor, bandpass=None, filter_name=None, + flux_to_counts_factor, exptimes=None, bandpass=None, filter_name=None, + wcs=None, rng=None, seed=None): """Add sources to an image. @@ -227,7 +228,7 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf, x & y positions of sources (pixel) at which sources should be added psf : galsim.Profile PSF for image - flux_to_counts_factor : float + flux_to_counts_factor : array_like physical fluxes in objlist (whether in profile SEDs or flux arrays) should be multiplied by this factor to convert to total counts in the image @@ -272,6 +273,10 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf, for i, obj in enumerate(objlist): t0 = time.time() image_pos = galsim.PositionD(xpos[i], ypos[i]) + if wcs is None: + pwcs = image.wcs.local(image_pos) + else: + pwcs = wcs.local(image_pos) profile = obj.profile if not chromatic: if obj.flux is None: @@ -282,21 +287,35 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf, psf0 = psf.at_position(xpos[i], ypos[i]) else: psf0 = psf - final = galsim.Convolve(profile * flux_to_counts_factor, psf0) + if isinstance(flux_to_counts_factor, list): + final = galsim.Convolve(profile * flux_to_counts_factor[i], psf0) + else: + final = galsim.Convolve(profile * flux_to_counts_factor, psf0) if chromatic: stamp = final.drawImage( - bandpass, center=image_pos, wcs=image.wcs.local(image_pos), + bandpass, center=image_pos, wcs=pwcs, method='phot', rng=rng) else: try: stamp = final.drawImage(center=image_pos, - wcs=image.wcs.local(image_pos)) + wcs=pwcs) except galsim.GalSimFFTSizeError: log.warning(f'Skipping source {i} due to too ' f'large FFT needed for desired accuracy.') - bounds = stamp.bounds & image.bounds + if exptimes is not None: + stamp /= exptimes[i] + + if isinstance(image, u.Quantity): + im_bounds = galsim.BoundsI(galsim.PositionI(0,0), galsim.PositionI(image.shape)) + else: + im_bounds = image.bounds + bounds = stamp.bounds & im_bounds if bounds.area() > 0: - image[bounds] += stamp[bounds] + if isinstance(image, u.Quantity): + image[bounds.getXMin():bounds.getXMax()+1, + bounds.getYMin():bounds.getYMax()+1] += stamp[bounds].array * image.unit + else: + image[bounds] += stamp[bounds] counts = np.sum(stamp[bounds].array) else: counts = 0 @@ -306,6 +325,105 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf, return outinfo +# def add_objects_to_image_old(image, objlist, xpos, ypos, psf, +# flux_to_counts_factor, bandpass=None, filter_name=None, +# rng=None, seed=None): +# """Add sources to an image. + +# Note: this includes Poisson noise when photon shooting is used +# (i.e., for chromatic source profiles), and otherwise is noise free. + +# Parameters +# ---------- +# image : galsim.Image +# Image to which to add sources with associated WCS. +# objlist : list[CatalogObject] +# Objects to add to image +# xpos, ypos : array_like +# x & y positions of sources (pixel) at which sources should be added +# psf : galsim.Profile +# PSF for image +# flux_to_counts_factor : float +# physical fluxes in objlist (whether in profile SEDs or flux arrays) +# should be multiplied by this factor to convert to total counts in the +# image +# bandpass : galsim.Bandpass +# bandpass in which image is being rendered. This is used only in cases +# where chromatic profiles & PSFs are being used. +# filter_name : str +# filter to use to select appropriate flux from objlist. This is only +# used when achromatic PSFs and sources are being rendered. +# rng : galsim.BaseDeviate +# random number generator to use +# seed : int +# seed to use for random number generator + +# Returns +# ------- +# outinfo : np.ndarray +# Array structure containing rows for each source. The columns give +# the total number of counts from the source entering the image and +# the time taken to render the source. +# """ +# if rng is None and seed is None: +# seed = 143 +# log.warning( +# 'No RNG set, constructing a new default RNG from default seed.') +# if rng is None: +# rng = galsim.UniformDeviate(seed) + +# log.info('Adding sources to image...') +# nrender = 0 + +# chromatic = False +# if len(objlist) > 0 and objlist[0].profile.spectral: +# chromatic = True +# if len(objlist) > 0 and chromatic and bandpass is None: +# raise ValueError('bandpass must be set for chromatic PSF rendering.') +# if len(objlist) > 0 and not chromatic and filter_name is None: +# raise ValueError('must specify filter when using achromatic PSF ' +# 'rendering.') + +# outinfo = np.zeros(len(objlist), dtype=[('counts', 'f4'), ('time', 'f4')]) +# for i, obj in enumerate(objlist): +# t0 = time.time() +# image_pos = galsim.PositionD(xpos[i], ypos[i]) +# profile = obj.profile +# if not chromatic: +# if obj.flux is None: +# raise ValueError('Non-chromatic sources must have specified ' +# 'fluxes!') +# profile = profile.withFlux(obj.flux[filter_name]) +# if hasattr(psf, 'at_position'): +# psf0 = psf.at_position(xpos[i], ypos[i]) +# else: +# psf0 = psf +# final = galsim.Convolve(profile * flux_to_counts_factor, psf0) +# if chromatic: +# stamp = final.drawImage( +# bandpass, center=image_pos, wcs=image.wcs.local(image_pos), +# method='phot', rng=rng) +# else: +# try: +# stamp = final.drawImage(center=image_pos, +# wcs=image.wcs.local(image_pos)) +# except galsim.GalSimFFTSizeError: +# log.warning(f'Skipping source {i} due to too ' +# f'large FFT needed for desired accuracy.') +# bounds = stamp.bounds & image.bounds +# if bounds.area() > 0: +# image[bounds] += stamp[bounds] +# counts = np.sum(stamp[bounds].array) +# else: +# counts = 0 +# nrender += 1 +# outinfo[i] = (counts, time.time() - t0) +# log.info('Rendered %d sources...' % nrender) +# return outinfo + + + + def simulate_counts_generic(image, exptime, objlist=None, psf=None, zpflux=None, sky=None, dark=None, diff --git a/romanisim/l3.py b/romanisim/l3.py index 84faaed2..de2c0cb5 100644 --- a/romanisim/l3.py +++ b/romanisim/l3.py @@ -26,7 +26,8 @@ CENTER_SCA = 2 -def add_objects_to_l3(l3_mos, xpos=None, ypos=None, coords=None, coords_unit='rad', wcs=None, psf=None, rng=None, seed=None): +def add_objects_to_l3(l3_mos, source_cat, xpos=None, ypos=None, coords=None, + coords_unit='rad', wcs=None, psf=None, rng=None, seed=None): """Add objects to a Level 3 mosaic Parameters @@ -53,18 +54,30 @@ def add_objects_to_l3(l3_mos, xpos=None, ypos=None, coords=None, coords_unit='ra # Create PSF psf = romanisim.psf.make_psf(filter_name=filter_name, sca=CENTER_SCA, chromatic=False, webbpsf=True) - if coords is None: + if any(pos is None for pos in [xpos,ypos]): + # Need to construct position arrays + if coords is None: + coords = np.array([[o.sky_pos.ra.rad, o.sky_pos.dec.rad] + for o in source_cat]) + coords_unit='rad' # Generate x,y positions for sources sourcecountsall = galsim.ImageF(l3_mos.data.shape[0], l3_mos.data.shape[1], wcs=wcs, xmin=0, ymin=0) xpos, ypos = sourcecountsall.wcs.radecToxy(coords[:, 0], coords[:, 1], coords_unit) - elif all(pos is not None for pos in [xpos,ypos]): - log.error('No (xpos and ypos) or coords set. Adding zero objects.') - return + + # if coords is not None: + # # Generate x,y positions for sources + # sourcecountsall = galsim.ImageF(l3_mos.data.shape[0], l3_mos.data.shape[1], wcs=wcs, xmin=0, ymin=0) + # xpos, ypos = sourcecountsall.wcs.radecToxy(coords[:, 0], coords[:, 1], coords_unit) + # elif all(pos is not None for pos in [xpos,ypos]): + # log.error('No (xpos and ypos) or coords set. Adding zero objects.') + # return + + print(f"XXX wcs = {wcs}") xpos_idx = [round(x) for x in xpos] ypos_idx = [round(y) for y in ypos] - if ((min(xpos_idx) < 0) or (min(ypos_idx)) or (max(xpos_idx) > l3_mos.data.shape[0]) + if ((min(xpos_idx) < 0) or (min(ypos_idx) < 0) or (max(xpos_idx) > l3_mos.data.shape[0]) or (max(ypos_idx) > l3_mos.data.shape[1])): log.error(f"A source is out of bounds! " f"Source min x,y = {min(xpos_idx)},{min(ypos_idx)}, max = {max(xpos_idx)},{max(ypos_idx)} " @@ -75,374 +88,46 @@ def add_objects_to_l3(l3_mos, xpos=None, ypos=None, coords=None, coords_unit='ra Ct_all = np.divide(l3_mos.data.value, l3_mos.var_poisson.value, out=np.ones(l3_mos.data.shape), where=l3_mos.var_poisson.value != 0) + Ct = [] # Cycle over sources and add them to the mosaic for idx, (x, y) in enumerate(zip(xpos_idx, ypos_idx)): # Set scaling factor for injected sources # Flux / sigma_p^2 if l3_mos.var_poisson[x][y].value != 0: - Ct = math.fabs(l3_mos.data[x][y].value / l3_mos.var_poisson[x][y].value) + Ct.append(math.fabs(l3_mos.data[x][y].value / l3_mos.var_poisson[x][y].value)) # elif l3_mos.data[x][y].value != 0: # Ct = math.fabs(l3_mos.data[x][y].value) else: - continue - - # Create empty postage stamp galsim source image - sourcecounts = galsim.ImageF(l3_mos.data.shape[0], l3_mos.data.shape[1], wcs=wcs, xmin=0, ymin=0) - - # Simulate source postage stamp - romanisim.image.add_objects_to_image(sourcecounts, [source_cat[idx]], xpos=[xpos[idx]], - ypos=[ypos[idx]], psf=psf, flux_to_counts_factor=Ct, - bandpass=[filter_name], filter_name=filter_name, rng=rng) - - # Scale the source image back by its flux ratios - sourcecounts /= Ct - - # Add sources to the original mosaic data array - l3_mos.data = (l3_mos.data.value + sourcecounts.array) * l3_mos.data.unit - # l3_mos.data = (l3_mos.data.value + np.swapaxes(sourcecounts.array, 0, 1)) * l3_mos.data.unit - - # Note for the future - other noise sources (read and flat) need to be implemented + Ct.append(1.0) + + # Create empty postage stamp galsim source image + # sourcecounts = galsim.ImageF(l3_mos.data.shape[0], l3_mos.data.shape[1], wcs=wcs, xmin=0, ymin=0) + # sourcecounts = galsim.ImageF(l3_mos.data.value, wcs=wcs, xmin=0, ymin=0) + # sourcecounts = galsim.ImageF(l3_mos.data, wcs=wcs, xmin=0, ymin=0) + + # unit = l3_mos.data.unit + + # Simulate source postage stamp + # romanisim.image.add_objects_to_image_old(sourcecounts, source_cat, xpos=xpos, + romanisim.image.add_objects_to_image(l3_mos.data, source_cat, xpos=xpos, + ypos=ypos, psf=psf, flux_to_counts_factor=Ct, exptimes=Ct, + bandpass=[filter_name], filter_name=filter_name, + wcs=wcs, rng=rng, seed=seed) + + # # Scale the source image back by its flux ratios + # sourcecounts /= Ct + + # Add sources to the original mosaic data array + # l3_mos.data = (l3_mos.data.value + sourcecounts) * l3_mos.data.unit + # l3_mos.data = (l3_mos.data.value + np.swapaxes(sourcecounts.array, 0, 1)) * l3_mos.data.unit + # print(f"XXX type(l3_mos.data.value) = {type(l3_mos.data.value)}") + # l3_mos.data = sourcecounts.array * unit + + + # Note for the future - other noise sources (read and flat) need to be implemented # Set new poisson variance - # l3_mos.var_poisson = l3_mos.data.value / Ct_all - l3_mos.var_poisson = np.divide(l3_mos.data.value, Ct_all, out=l3_mos.var_poisson.value, - where=(Ct_all != 0)) * l3_mos.var_poisson.unit + l3_mos.var_poisson = (l3_mos.data.value / Ct_all) * l3_mos.var_poisson.unit # l3_mos is updated in place, so no return return None - - -def make_l3(image, context, metadata, exptimes, rng=None, seed=None, - saturation=None, unit_factor=1): - """TBD - """ - - log.info('Apportioning counts to mosaic...') - # resultants, dq = apportion_counts_to_resultants( - # counts.array, tij, inv_linearity=inv_linearity, crparam=crparam, - # persistence=persistence, tstart=tstart, - # rng=rng, seed=seed) - - # roman.addReciprocityFailure(resultants_object) - - # code from apportion_counts_to_resultants - - # mosaic = np.swapaxes(image.array, 0, 1) - mosaic = image.array.copy() - - # Set rng for creating readnoise, flat noise - if rng is None and seed is None: - seed = 46 - log.warning( - 'No RNG set, constructing a new default RNG from default seed.') - if rng is None: - rng = galsim.GaussianDeviate(seed) - - # Need to add read noise, poisson noise, and flat noise - - # Simulate read and flat as random noise - # if read_noise is not None and not isinstance(read_noise, u.Quantity): - # read_noise = read_noise * u.DN - # log.warning('Making up units for read noise.') - # resultants are now in counts. - # read noise is in counts. - # log.info('Adding read noise...') - # resultants = add_read_noise_to_resultants( - # resultants, tij, rng=rng, seed=seed, - # read_noise=read_noise) - - # Improve this for proper 3D contexts - # ctx_2d = np.reshape(context[-2:], context.shape[-2:]) - # mosaic = np.divide(counts_arr, context[-2:], where=context[-2:]!=0, out=np.zeros(counts_arr.shape, dtype='f8')) - # mosaic = np.divide(counts_arr, ctx_2d, where=ctx_2d!=0, out=np.zeros(counts_arr.shape, dtype='f8')) - - # Mosaics may suffer saturation? If so rework the following. - # if saturation is None: - # saturation = parameters.reference_data['saturation'] - # mosaic = np.clip(mosaic, 0, 2 * 10**9).astype('i4') - # this maybe should be better applied at read time? - # it's not actually clear to me what the right thing to do - # is in detail. - # mosaic = np.clip(mosaic, 0 * u.DN, saturation) - # mosaic = np.clip(mosaic, 0, saturation.value) - # m = mosaic >= saturation.value - # dq = np.zeros(counts_arr.shape, dtype=np.uint32) - # dq[m] |= parameters.dqbits['saturated'] - - # Set mosaic to be a mosaic model - mosaic_mdl = maker_utils.mk_level3_mosaic(shape=mosaic.shape, meta=metadata) - - # Extract exposure time for each pixel - exptimes_pix = util.decode_context_times(context, exptimes) - - # Set data - mosaic_mdl.data = mosaic * unit_factor - mosaic_mdl.data = np.divide(mosaic_mdl.data.value, exptimes_pix, - out=np.zeros(mosaic_mdl.data.shape), - where=exptimes_pix != 0) * mosaic_mdl.data.unit - - # Context - # Binary index of images that contribute to the pixel - # Defined by geometry.. if not, all non zero = 1. - mosaic_mdl.context = context - - # Poisson noise - mosaic_mdl.var_poisson = unit_factor**2 / exptimes_pix**2 - - # Weight - # Use exptime weight - # --- - # DQ used if set above for saturation? - # bitvalue = interpret_bit_flags( - # bitvalue, flag_name_map={dq.name: dq.value for dq in pixel} ) - # if bitvalue is None: - # return np.ones(dqarr.shape, dtype=np.uint8) - # return np.logical_not(np.bitwise_and(dqarr, ~bitvalue)).astype(np.uint8) - # --- - # dqmask = build_mask(model.dq, good_bits) - # elif weight_type == "exptime": - # exptime = model.meta.exposure.exposure_time - # result = exptime * dqmask - - # Fill everything else: - # err, weight, var_rnoise, var_flat - - return mosaic_mdl - - -def generate_mosaic_geometry(): - """Generate a geometry map (context) for a mosaic dither pattern / set of pointings and roll angles - TBD - """ - pass - - -def simulate(metadata, cat, exptimes, context=None, - usecrds=True, webbpsf=True, seed=None, rng=None, - psf_keywords=dict(), **kwargs - ): - """TBD - """ - - 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.') - - # Create metadata object - meta = maker_utils.mk_mosaic_meta() - meta['wcs'] = None - - for key in parameters.default_mosaic_parameters_dictionary.keys(): - meta[key].update(parameters.default_mosaic_parameters_dictionary[key]) - - for key in metadata.keys(): - if key not in meta: - meta[key] = metadata[key] - # elif isinstance(meta[key], dict): - else: - meta[key].update(metadata[key]) - - # May need an equivalent of this this for L3 - # util.add_more_metadata(meta) - - filter_name = metadata['basic']['optical_element'] - - if rng is None and seed is None: - seed = 43 - log.warning( - 'No RNG set, constructing a new default RNG from default seed.') - if rng is None: - rng = galsim.UniformDeviate(seed) - - if context is None: - # Create geometry from the object list - twcs = romanisim.wcs.get_mosaic_wcs(meta) - - coords = np.array([[o.sky_pos.ra.rad, o.sky_pos.dec.rad] - for o in cat]) - - allx, ally = twcs.radecToxy(coords[:, 0], coords[:, 1], 'rad') - - # Obtain the sample extremums - xmin = min(allx) - xmax = max(allx) - ymin = min(ally) - ymax = max(ally) - - # Obtain WCS center - xcen, ycen = twcs.radecToxy(twcs.center.ra, twcs.center.dec, 'rad') - - # Determine maximum extremums from WCS center - xdiff = max([math.ceil(xmax - xcen), math.ceil(xcen - xmin)]) - ydiff = max([math.ceil(ymax - ycen), math.ceil(ycen - ymin)]) - - # Create context map preserving WCS center - context = np.ones((1, 2 * xdiff, 2 * ydiff), dtype=np.uint32) - - log.info('Simulating filter {0}...'.format(filter_name)) - mosaic, simcatobj = simulate_cps( - meta, cat, context, exptimes, rng=rng, usecrds=usecrds, - webbpsf=webbpsf, psf_keywords=psf_keywords) - - extras = {} - - # if reffiles: - # extras["simulate_reffiles"] = {} - # for key, value in reffiles.items(): - # extras["simulate_reffiles"][key] = value - - extras['simcatobj'] = simcatobj - - # TODO: Double check that this is a gwcs object - # extras['wcs'] = wcs.convert_wcs_to_gwcs(mosaic.wcs) - # extras['wcs'] = mosaic.wcs - - log.info('Simulation complete.') - # return mosaic, extras - return mosaic, extras - - -def simulate_cps(metadata, objlist, context, exptimes, - rng=None, seed=None, webbpsf=True, usecrds=False, - psf_keywords=dict()): - """TBD - """ - - if rng is None and seed is None: - seed = 43 - log.warning( - 'No RNG set, constructing a new default RNG from default seed.') - if rng is None: - rng = galsim.UniformDeviate(seed) - - filter_name = metadata['basic']['optical_element'] - - date = metadata['basic']['time_mean_mjd'] - - if not isinstance(date, astropy.time.Time): - date = astropy.time.Time(date, format='mjd') - - galsim_filter_name = romanisim.bandpass.roman2galsim_bandpass[filter_name] - bandpass = roman.getBandpasses(AB_zeropoint=True)[galsim_filter_name] - - # Generate WCS - moswcs = romanisim.wcs.get_mosaic_wcs(metadata, shape=context.shape[-2:]) - - # Is this needed? - chromatic = False - if (len(objlist) > 0 - and not isinstance(objlist, table.Table) # this case is always gray - and objlist[0].profile.spectral): - chromatic = True - - # With context, is possible that individual PSFs could be properly assigned? - # Would be a pain to implement - psf = romanisim.psf.make_psf(filter_name=filter_name, wcs=moswcs, sca=CENTER_SCA, - chromatic=chromatic, webbpsf=webbpsf, - variable=True, **psf_keywords) - - # Create initial galsim image - image = galsim.ImageF(context.shape[-2], context.shape[-1], wcs=moswcs, xmin=0, ymin=0) - - # Examine this in more detail to ensure it is correct - # Create base sky level - mos_cent_pos = moswcs.toWorld(image.true_center) - sky_level = roman.getSkyLevel(bandpass, world_pos=mos_cent_pos, - date=date.datetime, exptime=1) - sky_level *= (1.0 + roman.stray_light_fraction) - # sky_mosaic = image * 0 - # moswcs.makeSkyImage(sky_mosaic, sky_level) - # sky_mosaic += roman.thermal_backgrounds[galsim_filter_name] - abflux = romanisim.bandpass.get_abflux(filter_name) - - # Obtain physical unit conversion factor - unit_factor = parameters.reference_data['photom'][filter_name] - - # Maybe return a mosaic object as well? - mosaic, simcatobj = simulate_cps_generic( - image, metadata, context, exptimes, objlist=objlist, psf=psf, zpflux=abflux, sky=sky_level, - wcs=moswcs, rng=rng, seed=seed, unit_factor=unit_factor) - - return mosaic, simcatobj - - -def simulate_cps_generic(image, metadata, context, exptimes, objlist=None, psf=None, - zpflux=None, wcs=None, xpos=None, ypos=None, sky=None, - flat=None, rng=None, seed=None, unit_factor=1): - """TBD - """ - if rng is None and seed is None: - seed = 144 - log.warning( - 'No RNG set, constructing a new default RNG from default seed.') - if rng is None: - rng = galsim.UniformDeviate(seed) - - # Objlist wcs check - if (objlist is not None and len(objlist) > 0 - and image.wcs is None and (xpos is None or ypos is None)): - raise ValueError('xpos and ypos must be set if rendering objects ' - 'without a WCS.') - if objlist is None: - objlist = [] - if len(objlist) > 0 and xpos is None: - if isinstance(objlist, table.Table): - xpos, ypos = image.wcs.radecToxy( - np.radians(objlist['ra']), np.radians(objlist['dec']), 'rad') - else: - coord = np.array([[o.sky_pos.ra.rad, o.sky_pos.dec.rad] - for o in objlist]) - xpos, ypos = image.wcs.radecToxy(coord[:, 0], coord[:, 1], 'rad') - # use private vectorized transformation - if xpos is not None: - xpos = np.array(xpos) - if ypos is not None: - ypos = np.array(ypos) - - if len(objlist) > 0 and psf is None: - raise ValueError('Must provide a PSF if you want to render objects.') - - if flat is None: - flat = 1.0 - - # # for some reason, galsim doesn't like multiplying an SED by 1, but it's - # # okay with multiplying an SED by 1.0, so we cast to float. - 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. - - chromatic = False - if len(objlist) > 0 and objlist[0].profile.spectral: - chromatic = True - - # add Poisson noise if we made a noiseless, not-photon-shooting - # image. - poisson_noise = galsim.PoissonNoise(rng, sky_level=sky) - if not chromatic: - image.addNoise(poisson_noise) - - # Add Objects requires a mosaic image - # Maybe violates the notion of a "generic" count simulator? - mosaic_mdl = make_l3(image, context, metadata, exptimes, unit_factor=unit_factor) - - if wcs is not None: - mosaic_mdl['wcs'] = wcs - - # Add objects to mosaic model - objinfo = add_objects_to_l3( - mosaic_mdl, objlist, rng=rng) - objinfo = np.zeros( - len(objlist), - dtype=[('x', 'f4'), ('y', 'f4'), ('counts', 'f4'), ('time', 'f4')]) - if len(objlist) > 0: - objinfo['x'] = xpos - objinfo['y'] = ypos - - return mosaic_mdl, objinfo diff --git a/romanisim/tests/test_l3.py b/romanisim/tests/test_l3.py index 6cc580ef..db8f6669 100644 --- a/romanisim/tests/test_l3.py +++ b/romanisim/tests/test_l3.py @@ -85,12 +85,29 @@ def test_inject_sources_into_mosaic(): sc_table['ra'][3], sc_table['dec'][3] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value source_cat = catalog.table_to_catalog(sc_table, ["F158"]) + coords = np.array([[o.sky_pos.ra.rad, o.sky_pos.dec.rad] + for o in source_cat]) # Copy original Mosaic before adding sources l3_mos_orig = l3_mos.copy() + l3_mos_orig.data = l3_mos.data.copy() + l3_mos_orig.var_poisson = l3_mos.var_poisson.copy() # Add source_cat objects to mosaic l3.add_objects_to_l3(l3_mos, source_cat, seed=rng_seed) + # l3.add_objects_to_l3(l3_mos, source_cat, coords=coords, seed=rng_seed) + + import plotly.express as px + + fig1 = px.imshow(l3_mos_orig.data.value, title='Orig Mosaic Data', labels={'color': 'MJy / sr'}) + fig1.show() + + fig2 = px.imshow(l3_mos.data.value, title='Injected Mosaic Data', labels={'color': 'MJy / sr'}) + fig2.show() + + fig3 = px.imshow((l3_mos.data.value - l3_mos_orig.data.value), title='Diff Mosaic Data', labels={'color': 'MJy / sr'}) + fig3.show() + # Ensure that every data pixel value has increased or # remained the same with the new sources injected @@ -100,6 +117,11 @@ def test_inject_sources_into_mosaic(): # remained the same with the new sources injected # Numpy isclose is needed to determine equality, due to float precision issues close_mask = np.isclose(l3_mos.var_poisson.value, l3_mos_orig.var_poisson.value, rtol=1e-06) + + fig4 = px.imshow(close_mask.astype(int), title='Makes', labels={'color': 'T/F'}) + fig4.show() + + assert False in close_mask assert np.all(l3_mos.var_poisson.value[~close_mask] > l3_mos_orig.var_poisson.value[~close_mask]) # Create log entry and artifacts @@ -113,48 +135,3 @@ def test_inject_sources_into_mosaic(): 'source_cat_table': sc_table, } af.write_to(os.path.join(artifactdir, 'dms232.asdf')) - - -def test_sim_mosaic(): - """Simulating mosaic from catalog file. - """ - - ra_ref = 1.0 - dec_ref = 4.0 - - metadata = copy.deepcopy(parameters.default_mosaic_parameters_dictionary) - metadata['basic']['optical_element'] = 'F158' - metadata['wcsinfo']['ra_ref'] = ra_ref - metadata['wcsinfo']['dec_ref'] = dec_ref - metadata['wcsinfo']['pixel_scale'] = 0.11 - metadata['wcsinfo']['pixel_scale_local'] = 0.11 - metadata['wcsinfo']['v2_ref'] = 0 - metadata['wcsinfo']['v3_ref'] = 0 - - exptimes = [600] - - cen = SkyCoord(ra=ra_ref * u.deg, dec=dec_ref * u.deg) - cat = catalog.make_dummy_table_catalog(cen, radius=0.01, nobj=100) - cat['F158'] = cat['F158'] * 10e10 - - source_cat = catalog.table_to_catalog(cat, ["F158"]) - - mosaic, extras = l3.simulate(metadata, source_cat, exptimes) - - import plotly.express as px - - fig1 = px.imshow(mosaic.data.value, title='Mosaic Data', labels={'color': 'MJy / sr'}) - fig1.show() - - fig2 = px.imshow(mosaic.context[-2:].reshape(mosaic.context.shape[-2:]), title='Mosaic Context', labels={'color': 'File Number'}) - fig2.show() - - pos_vals = mosaic.data.value.copy() - pos_vals[pos_vals <= 0] = 0.00000000001 - - fig3 = px.imshow(np.log(pos_vals), title='Mosaic Data (log)', labels={'color': 'MJy / sr'}) - fig3.show() - -# TBD: Test with more complex context - -# TBD: Test of geometry construction From 6536f8e7b8dc8cd24055da79edec5f5d81a89bf9 Mon Sep 17 00:00:00 2001 From: PaulHuwe Date: Sun, 7 Jul 2024 17:17:23 -0400 Subject: [PATCH 5/9] Initial PR --- romanisim/image.py | 121 ++++--------------------------------- romanisim/l3.py | 98 +++++++++--------------------- romanisim/tests/test_l3.py | 50 +++++++-------- 3 files changed, 64 insertions(+), 205 deletions(-) diff --git a/romanisim/image.py b/romanisim/image.py index d70c03cd..074d3cc8 100644 --- a/romanisim/image.py +++ b/romanisim/image.py @@ -210,9 +210,9 @@ def trim_objlist(objlist, image): def add_objects_to_image(image, objlist, xpos, ypos, psf, - flux_to_counts_factor, exptimes=None, bandpass=None, filter_name=None, - wcs=None, - rng=None, seed=None): + flux_to_counts_factor, exptimes=None, + bandpass=None, filter_name=None, + wcs=None, rng=None, seed=None): """Add sources to an image. Note: this includes Poisson noise when photon shooting is used @@ -220,24 +220,28 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf, Parameters ---------- - image : galsim.Image - Image to which to add sources with associated WCS. + image : galsim.Image or astropy.Quantity[float] + Image to which to add sources with associated WCS. Updated in place. objlist : list[CatalogObject] Objects to add to image xpos, ypos : array_like x & y positions of sources (pixel) at which sources should be added psf : galsim.Profile PSF for image - flux_to_counts_factor : array_like + flux_to_counts_factor : float or list physical fluxes in objlist (whether in profile SEDs or flux arrays) should be multiplied by this factor to convert to total counts in the image + exptimes: array_like + Exposure times to scale back to rate units bandpass : galsim.Bandpass bandpass in which image is being rendered. This is used only in cases where chromatic profiles & PSFs are being used. filter_name : str filter to use to select appropriate flux from objlist. This is only used when achromatic PSFs and sources are being rendered. + wcs : galsim.GSFitsWCS + WCS corresponding to image rng : galsim.BaseDeviate random number generator to use seed : int @@ -306,14 +310,14 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf, stamp /= exptimes[i] if isinstance(image, u.Quantity): - im_bounds = galsim.BoundsI(galsim.PositionI(0,0), galsim.PositionI(image.shape)) + im_bounds = galsim.BoundsI(galsim.PositionI(0, 0), galsim.PositionI(image.shape)) else: im_bounds = image.bounds bounds = stamp.bounds & im_bounds if bounds.area() > 0: if isinstance(image, u.Quantity): - image[bounds.getXMin():bounds.getXMax()+1, - bounds.getYMin():bounds.getYMax()+1] += stamp[bounds].array * image.unit + image[bounds.getXMin():(bounds.getXMax() + 1), + bounds.getYMin():(bounds.getYMax() + 1)] += stamp[bounds].array * image.unit else: image[bounds] += stamp[bounds] counts = np.sum(stamp[bounds].array) @@ -325,105 +329,6 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf, return outinfo -# def add_objects_to_image_old(image, objlist, xpos, ypos, psf, -# flux_to_counts_factor, bandpass=None, filter_name=None, -# rng=None, seed=None): -# """Add sources to an image. - -# Note: this includes Poisson noise when photon shooting is used -# (i.e., for chromatic source profiles), and otherwise is noise free. - -# Parameters -# ---------- -# image : galsim.Image -# Image to which to add sources with associated WCS. -# objlist : list[CatalogObject] -# Objects to add to image -# xpos, ypos : array_like -# x & y positions of sources (pixel) at which sources should be added -# psf : galsim.Profile -# PSF for image -# flux_to_counts_factor : float -# physical fluxes in objlist (whether in profile SEDs or flux arrays) -# should be multiplied by this factor to convert to total counts in the -# image -# bandpass : galsim.Bandpass -# bandpass in which image is being rendered. This is used only in cases -# where chromatic profiles & PSFs are being used. -# filter_name : str -# filter to use to select appropriate flux from objlist. This is only -# used when achromatic PSFs and sources are being rendered. -# rng : galsim.BaseDeviate -# random number generator to use -# seed : int -# seed to use for random number generator - -# Returns -# ------- -# outinfo : np.ndarray -# Array structure containing rows for each source. The columns give -# the total number of counts from the source entering the image and -# the time taken to render the source. -# """ -# if rng is None and seed is None: -# seed = 143 -# log.warning( -# 'No RNG set, constructing a new default RNG from default seed.') -# if rng is None: -# rng = galsim.UniformDeviate(seed) - -# log.info('Adding sources to image...') -# nrender = 0 - -# chromatic = False -# if len(objlist) > 0 and objlist[0].profile.spectral: -# chromatic = True -# if len(objlist) > 0 and chromatic and bandpass is None: -# raise ValueError('bandpass must be set for chromatic PSF rendering.') -# if len(objlist) > 0 and not chromatic and filter_name is None: -# raise ValueError('must specify filter when using achromatic PSF ' -# 'rendering.') - -# outinfo = np.zeros(len(objlist), dtype=[('counts', 'f4'), ('time', 'f4')]) -# for i, obj in enumerate(objlist): -# t0 = time.time() -# image_pos = galsim.PositionD(xpos[i], ypos[i]) -# profile = obj.profile -# if not chromatic: -# if obj.flux is None: -# raise ValueError('Non-chromatic sources must have specified ' -# 'fluxes!') -# profile = profile.withFlux(obj.flux[filter_name]) -# if hasattr(psf, 'at_position'): -# psf0 = psf.at_position(xpos[i], ypos[i]) -# else: -# psf0 = psf -# final = galsim.Convolve(profile * flux_to_counts_factor, psf0) -# if chromatic: -# stamp = final.drawImage( -# bandpass, center=image_pos, wcs=image.wcs.local(image_pos), -# method='phot', rng=rng) -# else: -# try: -# stamp = final.drawImage(center=image_pos, -# wcs=image.wcs.local(image_pos)) -# except galsim.GalSimFFTSizeError: -# log.warning(f'Skipping source {i} due to too ' -# f'large FFT needed for desired accuracy.') -# bounds = stamp.bounds & image.bounds -# if bounds.area() > 0: -# image[bounds] += stamp[bounds] -# counts = np.sum(stamp[bounds].array) -# else: -# counts = 0 -# nrender += 1 -# outinfo[i] = (counts, time.time() - t0) -# log.info('Rendered %d sources...' % nrender) -# return outinfo - - - - def simulate_counts_generic(image, exptime, objlist=None, psf=None, zpflux=None, sky=None, dark=None, diff --git a/romanisim/l3.py b/romanisim/l3.py index de2c0cb5..d4653745 100644 --- a/romanisim/l3.py +++ b/romanisim/l3.py @@ -26,7 +26,7 @@ CENTER_SCA = 2 -def add_objects_to_l3(l3_mos, source_cat, xpos=None, ypos=None, coords=None, +def add_objects_to_l3(l3_mos, source_cat, exptimes, xpos=None, ypos=None, coords=None, coords_unit='rad', wcs=None, psf=None, rng=None, seed=None): """Add objects to a Level 3 mosaic @@ -36,6 +36,22 @@ def add_objects_to_l3(l3_mos, source_cat, xpos=None, ypos=None, coords=None, Mosaic of images source_cat : list List of catalog objects to add to l3_mos + exptimes : list + Exposure times to scale back to rate units + xpos, ypos : array_like + x & y positions of sources (pixel) at which sources should be added + coords : array_like + ra & dec positions of sources (coords_unit) at which sources should be added + coords_unit : string + units of coords + wcs : galsim.GSFitsWCS + WCS corresponding to image + psf : galsim.Profile + PSF for image + rng : galsim.BaseDeviate + random number generator to use + seed : int + seed to use for random number generator Returns ------- @@ -46,88 +62,30 @@ def add_objects_to_l3(l3_mos, source_cat, xpos=None, ypos=None, coords=None, # Obtain optical element filter_name = l3_mos.meta.basic.optical_element + # Generate WCS (if needed) if wcs is None: - # Generate WCS wcs = romanisim.wcs.get_mosaic_wcs(l3_mos.meta, shape=l3_mos.data.shape) + # Create PSF (if needed) if psf is None: - # Create PSF psf = romanisim.psf.make_psf(filter_name=filter_name, sca=CENTER_SCA, chromatic=False, webbpsf=True) - if any(pos is None for pos in [xpos,ypos]): - # Need to construct position arrays + # Create position arrays (if needed) + if any(pos is None for pos in [xpos, ypos]): + # Create coordinates (if needed) if coords is None: coords = np.array([[o.sky_pos.ra.rad, o.sky_pos.dec.rad] - for o in source_cat]) - coords_unit='rad' + for o in source_cat]) + coords_unit = 'rad' # Generate x,y positions for sources sourcecountsall = galsim.ImageF(l3_mos.data.shape[0], l3_mos.data.shape[1], wcs=wcs, xmin=0, ymin=0) xpos, ypos = sourcecountsall.wcs.radecToxy(coords[:, 0], coords[:, 1], coords_unit) - # if coords is not None: - # # Generate x,y positions for sources - # sourcecountsall = galsim.ImageF(l3_mos.data.shape[0], l3_mos.data.shape[1], wcs=wcs, xmin=0, ymin=0) - # xpos, ypos = sourcecountsall.wcs.radecToxy(coords[:, 0], coords[:, 1], coords_unit) - # elif all(pos is not None for pos in [xpos,ypos]): - # log.error('No (xpos and ypos) or coords set. Adding zero objects.') - # return - - print(f"XXX wcs = {wcs}") - - xpos_idx = [round(x) for x in xpos] - ypos_idx = [round(y) for y in ypos] - - if ((min(xpos_idx) < 0) or (min(ypos_idx) < 0) or (max(xpos_idx) > l3_mos.data.shape[0]) - or (max(ypos_idx) > l3_mos.data.shape[1])): - log.error(f"A source is out of bounds! " - f"Source min x,y = {min(xpos_idx)},{min(ypos_idx)}, max = {max(xpos_idx)},{max(ypos_idx)} " - f"Image min x,y = {0},{0}, max = {l3_mos.data.shape[0]},{l3_mos.data.shape[1]}") - - # Create overall scaling factor map - # Ct_all = (l3_mos.data.value / l3_mos.var_poisson) - Ct_all = np.divide(l3_mos.data.value, l3_mos.var_poisson.value, - out=np.ones(l3_mos.data.shape), where=l3_mos.var_poisson.value != 0) - - Ct = [] - # Cycle over sources and add them to the mosaic - for idx, (x, y) in enumerate(zip(xpos_idx, ypos_idx)): - # Set scaling factor for injected sources - # Flux / sigma_p^2 - if l3_mos.var_poisson[x][y].value != 0: - Ct.append(math.fabs(l3_mos.data[x][y].value / l3_mos.var_poisson[x][y].value)) - # elif l3_mos.data[x][y].value != 0: - # Ct = math.fabs(l3_mos.data[x][y].value) - else: - Ct.append(1.0) - - # Create empty postage stamp galsim source image - # sourcecounts = galsim.ImageF(l3_mos.data.shape[0], l3_mos.data.shape[1], wcs=wcs, xmin=0, ymin=0) - # sourcecounts = galsim.ImageF(l3_mos.data.value, wcs=wcs, xmin=0, ymin=0) - # sourcecounts = galsim.ImageF(l3_mos.data, wcs=wcs, xmin=0, ymin=0) - - # unit = l3_mos.data.unit - - # Simulate source postage stamp - # romanisim.image.add_objects_to_image_old(sourcecounts, source_cat, xpos=xpos, - romanisim.image.add_objects_to_image(l3_mos.data, source_cat, xpos=xpos, - ypos=ypos, psf=psf, flux_to_counts_factor=Ct, exptimes=Ct, - bandpass=[filter_name], filter_name=filter_name, - wcs=wcs, rng=rng, seed=seed) - - # # Scale the source image back by its flux ratios - # sourcecounts /= Ct - # Add sources to the original mosaic data array - # l3_mos.data = (l3_mos.data.value + sourcecounts) * l3_mos.data.unit - # l3_mos.data = (l3_mos.data.value + np.swapaxes(sourcecounts.array, 0, 1)) * l3_mos.data.unit - # print(f"XXX type(l3_mos.data.value) = {type(l3_mos.data.value)}") - # l3_mos.data = sourcecounts.array * unit - - - # Note for the future - other noise sources (read and flat) need to be implemented - - # Set new poisson variance - l3_mos.var_poisson = (l3_mos.data.value / Ct_all) * l3_mos.var_poisson.unit + romanisim.image.add_objects_to_image(l3_mos.data, source_cat, xpos=xpos, ypos=ypos, + psf=psf, flux_to_counts_factor=exptimes, exptimes=exptimes, + bandpass=[filter_name], filter_name=filter_name, + wcs=wcs, rng=rng, seed=seed) # l3_mos is updated in place, so no return return None diff --git a/romanisim/tests/test_l3.py b/romanisim/tests/test_l3.py index db8f6669..669b2860 100644 --- a/romanisim/tests/test_l3.py +++ b/romanisim/tests/test_l3.py @@ -4,6 +4,7 @@ import os import copy +import math import numpy as np import galsim from galsim import roman @@ -75,39 +76,38 @@ def test_inject_sources_into_mosaic(): "half_light_radius": 4 * [0.0], "pa": 4 * [0.0], "ba": 4 * [1.0], "F158": 4 * [1.0]} sc_table = table.Table(sc_dict) - xpos, ypos = 50, 50 - sc_table["ra"][0], sc_table["dec"][0] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value - xpos, ypos = 50, 150 - sc_table['ra'][1], sc_table['dec'][1] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value - xpos, ypos = 150, 50 - sc_table['ra'][2], sc_table['dec'][2] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value - xpos, ypos = 150, 150 - sc_table['ra'][3], sc_table['dec'][3] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value + # Set locations + xpos_idx = [50, 50, 150, 150] + ypos_idx = [50, 150, 50, 150] + + # Populate flux scaling ratio and catalog + Ct = [] + for idx, (x, y) in enumerate(zip(xpos_idx, ypos_idx)): + # Set scaling factor for injected sources + # Flux / sigma_p^2 + if l3_mos.var_poisson[x][y].value != 0: + Ct.append(math.fabs(l3_mos.data[x][y].value / l3_mos.var_poisson[x][y].value)) + else: + Ct.append(1.0) + + sc_table["ra"][idx], sc_table["dec"][idx] = (twcs._radec(x, y) * u.rad).to(u.deg).value source_cat = catalog.table_to_catalog(sc_table, ["F158"]) - coords = np.array([[o.sky_pos.ra.rad, o.sky_pos.dec.rad] - for o in source_cat]) - # Copy original Mosaic before adding sources + # Copy original Mosaic before adding sources as sources are added in place l3_mos_orig = l3_mos.copy() l3_mos_orig.data = l3_mos.data.copy() l3_mos_orig.var_poisson = l3_mos.var_poisson.copy() # Add source_cat objects to mosaic - l3.add_objects_to_l3(l3_mos, source_cat, seed=rng_seed) - # l3.add_objects_to_l3(l3_mos, source_cat, coords=coords, seed=rng_seed) - - import plotly.express as px - - fig1 = px.imshow(l3_mos_orig.data.value, title='Orig Mosaic Data', labels={'color': 'MJy / sr'}) - fig1.show() + l3.add_objects_to_l3(l3_mos, source_cat, Ct, seed=rng_seed) - fig2 = px.imshow(l3_mos.data.value, title='Injected Mosaic Data', labels={'color': 'MJy / sr'}) - fig2.show() - - fig3 = px.imshow((l3_mos.data.value - l3_mos_orig.data.value), title='Diff Mosaic Data', labels={'color': 'MJy / sr'}) - fig3.show() + # Create overall scaling factor map + Ct_all = np.divide(l3_mos_orig.data.value, l3_mos_orig.var_poisson.value, + out=np.ones(l3_mos_orig.data.shape), where=l3_mos_orig.var_poisson.value != 0) + # Set new poisson variance + l3_mos.var_poisson = (l3_mos.data.value / Ct_all) * l3_mos.var_poisson.unit # Ensure that every data pixel value has increased or # remained the same with the new sources injected @@ -117,10 +117,6 @@ def test_inject_sources_into_mosaic(): # remained the same with the new sources injected # Numpy isclose is needed to determine equality, due to float precision issues close_mask = np.isclose(l3_mos.var_poisson.value, l3_mos_orig.var_poisson.value, rtol=1e-06) - - fig4 = px.imshow(close_mask.astype(int), title='Makes', labels={'color': 'T/F'}) - fig4.show() - assert False in close_mask assert np.all(l3_mos.var_poisson.value[~close_mask] > l3_mos_orig.var_poisson.value[~close_mask]) From 65898b9a5b2e32f19e0858c5911c1bd3fd70106c Mon Sep 17 00:00:00 2001 From: PaulHuwe Date: Wed, 10 Jul 2024 11:40:01 -0400 Subject: [PATCH 6/9] Addressed PR comments and added flake8 fixes. --- romanisim/image.py | 14 +++----------- romanisim/l3.py | 17 ++++++++++++----- romanisim/psf.py | 8 ++++---- romanisim/tests/test_image.py | 3 +-- romanisim/tests/test_l1.py | 2 +- romanisim/tests/test_l3.py | 11 ++++++++--- romanisim/wcs.py | 8 ++++---- setup.py | 2 +- 8 files changed, 34 insertions(+), 31 deletions(-) diff --git a/romanisim/image.py b/romanisim/image.py index 3251e26c..9bc062e1 100644 --- a/romanisim/image.py +++ b/romanisim/image.py @@ -224,7 +224,7 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf, Parameters ---------- - image : galsim.Image or astropy.Quantity[float] + image : galsim.Image Image to which to add sources with associated WCS. Updated in place. objlist : list[CatalogObject] Objects to add to image @@ -313,17 +313,9 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf, if exptimes is not None: stamp /= exptimes[i] - if isinstance(image, u.Quantity): - im_bounds = galsim.BoundsI(galsim.PositionI(0, 0), galsim.PositionI(image.shape)) - else: - im_bounds = image.bounds - bounds = stamp.bounds & im_bounds + bounds = stamp.bounds & image.bounds if bounds.area() > 0: - if isinstance(image, u.Quantity): - image[bounds.getXMin():(bounds.getXMax() + 1), - bounds.getYMin():(bounds.getYMax() + 1)] += stamp[bounds].array * image.unit - else: - image[bounds] += stamp[bounds] + image[bounds] += stamp[bounds] counts = np.sum(stamp[bounds].array) else: counts = 0 diff --git a/romanisim/l3.py b/romanisim/l3.py index d4653745..f6f0bbd8 100644 --- a/romanisim/l3.py +++ b/romanisim/l3.py @@ -26,7 +26,7 @@ CENTER_SCA = 2 -def add_objects_to_l3(l3_mos, source_cat, exptimes, xpos=None, ypos=None, coords=None, +def add_objects_to_l3(l3_mos, source_cat, exptimes, xpos=None, ypos=None, coords=None, unit_factor=1.0, coords_unit='rad', wcs=None, psf=None, rng=None, seed=None): """Add objects to a Level 3 mosaic @@ -42,6 +42,8 @@ def add_objects_to_l3(l3_mos, source_cat, exptimes, xpos=None, ypos=None, coords x & y positions of sources (pixel) at which sources should be added coords : array_like ra & dec positions of sources (coords_unit) at which sources should be added + unit_factor: float + Factor to convert data to MJy / sr coords_unit : string units of coords wcs : galsim.GSFitsWCS @@ -70,6 +72,9 @@ def add_objects_to_l3(l3_mos, source_cat, exptimes, xpos=None, ypos=None, coords if psf is None: psf = romanisim.psf.make_psf(filter_name=filter_name, sca=CENTER_SCA, chromatic=False, webbpsf=True) + # Create Image canvas to add objects to + sourcecountsall = galsim.ImageF(l3_mos.data.value, wcs=wcs, xmin=0, ymin=0) + # Create position arrays (if needed) if any(pos is None for pos in [xpos, ypos]): # Create coordinates (if needed) @@ -78,14 +83,16 @@ def add_objects_to_l3(l3_mos, source_cat, exptimes, xpos=None, ypos=None, coords for o in source_cat]) coords_unit = 'rad' # Generate x,y positions for sources - sourcecountsall = galsim.ImageF(l3_mos.data.shape[0], l3_mos.data.shape[1], wcs=wcs, xmin=0, ymin=0) xpos, ypos = sourcecountsall.wcs.radecToxy(coords[:, 0], coords[:, 1], coords_unit) # Add sources to the original mosaic data array - romanisim.image.add_objects_to_image(l3_mos.data, source_cat, xpos=xpos, ypos=ypos, - psf=psf, flux_to_counts_factor=exptimes, exptimes=exptimes, - bandpass=[filter_name], filter_name=filter_name, + romanisim.image.add_objects_to_image(sourcecountsall, source_cat, xpos=xpos, ypos=ypos, + psf=psf, flux_to_counts_factor=[xpt * unit_factor for xpt in exptimes], + exptimes=exptimes, bandpass=[filter_name], filter_name=filter_name, wcs=wcs, rng=rng, seed=seed) + # Save array with added sources + l3_mos.data = sourcecountsall.array * l3_mos.data.unit + # l3_mos is updated in place, so no return return None diff --git a/romanisim/psf.py b/romanisim/psf.py index 40294500..a669751f 100644 --- a/romanisim/psf.py +++ b/romanisim/psf.py @@ -201,8 +201,8 @@ def at_position(self, x, y): # x = [0, off] -> 1 # x = [npix, infinity] -> 0 # linearly between those, likewise for y. - out = (self.psf['ll'] * wleft * wlow + - self.psf['lr'] * (1 - wleft) * wlow + - self.psf['ul'] * wleft * (1 - wlow) + - self.psf['ur'] * (1 - wleft) * (1 - wlow)) + out = (self.psf['ll'] * wleft * wlow + + self.psf['lr'] * (1 - wleft) * wlow + + self.psf['ul'] * wleft * (1 - wlow) + + self.psf['ur'] * (1 - wleft) * (1 - wlow)) return out diff --git a/romanisim/tests/test_image.py b/romanisim/tests/test_image.py index c56378f2..99319e85 100644 --- a/romanisim/tests/test_image.py +++ b/romanisim/tests/test_image.py @@ -18,7 +18,7 @@ import numpy as np import galsim from galsim import roman -from romanisim import image, parameters, catalog, psf, util, wcs, persistence, ramp, l1, l3 +from romanisim import image, parameters, catalog, psf, util, wcs, persistence, ramp, l1 from astropy.coordinates import SkyCoord from astropy import units as u from astropy.time import Time @@ -30,7 +30,6 @@ from metrics_logger.decorators import metrics_logger from romanisim import log from roman_datamodels.stnode import WfiScienceRaw, WfiImage -import roman_datamodels.maker_utils as maker_utils import romanisim.bandpass diff --git a/romanisim/tests/test_l1.py b/romanisim/tests/test_l1.py index d4a7aea4..5f944872 100644 --- a/romanisim/tests/test_l1.py +++ b/romanisim/tests/test_l1.py @@ -29,7 +29,7 @@ read_pattern_list = [ [[1 + x for x in range(10)], [11], [12 + x for x in range(10)], [30], - [40 + x for x in range(5)], [50 + x for x in range(100)]], + [40 + x for x in range(5)], [50 + x for x in range(100)]], [[1]], [[1], [10]], ] diff --git a/romanisim/tests/test_l3.py b/romanisim/tests/test_l3.py index 669b2860..34a7838c 100644 --- a/romanisim/tests/test_l3.py +++ b/romanisim/tests/test_l3.py @@ -59,11 +59,15 @@ def test_inject_sources_into_mosaic(): if key in l3_mos.meta: l3_mos.meta[key].update(metadata[key]) + # Obtain unit conversion factor + unit_factor = parameters.reference_data['photom'][l3_mos.meta.basic.optical_element] + # Populate the mosaic data array with gaussian noise from generators g1.generate(l3_mos.data.value[0:100, 0:100]) g2.generate(l3_mos.data.value[0:100, 100:200]) g3.generate(l3_mos.data.value[100:200, 0:100]) g4.generate(l3_mos.data.value[100:200, 100:200]) + l3_mos.data *= unit_factor.value # Define Poisson Noise of mosaic l3_mos.var_poisson.value[0:100, 0:100] = 0.01**2 @@ -85,8 +89,8 @@ def test_inject_sources_into_mosaic(): for idx, (x, y) in enumerate(zip(xpos_idx, ypos_idx)): # Set scaling factor for injected sources # Flux / sigma_p^2 - if l3_mos.var_poisson[x][y].value != 0: - Ct.append(math.fabs(l3_mos.data[x][y].value / l3_mos.var_poisson[x][y].value)) + if l3_mos.var_poisson[y, x].value != 0: + Ct.append(math.fabs(l3_mos.data[y, x].value / l3_mos.var_poisson[y, x].value)) else: Ct.append(1.0) @@ -100,7 +104,7 @@ def test_inject_sources_into_mosaic(): l3_mos_orig.var_poisson = l3_mos.var_poisson.copy() # Add source_cat objects to mosaic - l3.add_objects_to_l3(l3_mos, source_cat, Ct, seed=rng_seed) + l3.add_objects_to_l3(l3_mos, source_cat, Ct, unit_factor=unit_factor.value, seed=rng_seed) # Create overall scaling factor map Ct_all = np.divide(l3_mos_orig.data.value, l3_mos_orig.var_poisson.value, @@ -117,6 +121,7 @@ def test_inject_sources_into_mosaic(): # remained the same with the new sources injected # Numpy isclose is needed to determine equality, due to float precision issues close_mask = np.isclose(l3_mos.var_poisson.value, l3_mos_orig.var_poisson.value, rtol=1e-06) + assert False in close_mask assert np.all(l3_mos.var_poisson.value[~close_mask] > l3_mos_orig.var_poisson.value[~close_mask]) diff --git a/romanisim/wcs.py b/romanisim/wcs.py index 06d6a993..59c31380 100644 --- a/romanisim/wcs.py +++ b/romanisim/wcs.py @@ -220,15 +220,15 @@ def make_wcs(targ_pos, # V2V3 are in arcseconds, while SphericalToCartesian expects degrees, # so again start by scaling by 3600 tel2sky = ((models.Scale(1 / 3600) & models.Scale(1 / 3600)) - | gwcs.geometry.SphericalToCartesian(wrap_lon_at=wrap_v2_at) - | rot - | gwcs.geometry.CartesianToSpherical(wrap_lon_at=wrap_lon_at)) + | gwcs.geometry.SphericalToCartesian(wrap_lon_at=wrap_v2_at) + | rot + | gwcs.geometry.CartesianToSpherical(wrap_lon_at=wrap_lon_at)) tel2sky.name = 'v23tosky' detector = cf.Frame2D(name='detector', axes_order=(0, 1), unit=(u.pix, u.pix)) v2v3 = cf.Frame2D(name="v2v3", axes_order=(0, 1), - axes_names=("v2", "v3"), unit=(u.arcsec, u.arcsec)) + axes_names=("v2", "v3"), unit=(u.arcsec, u.arcsec)) world = cf.CelestialFrame(reference_frame=astropy.coordinates.ICRS(), name='world') diff --git a/setup.py b/setup.py index d2bb86d4..b9c44806 100644 --- a/setup.py +++ b/setup.py @@ -8,7 +8,7 @@ Options.docstrings = True Options.annotate = False -# importing these extension modules is tested in `.github/workflows/build.yml`; +# importing these extension modules is tested in `.github/workflows/build.yml`; # when adding new modules here, make sure to add them to the `test_command` entry there extensions = [Extension('romanisim.ramp_fit_casertano', ['romanisim/ramp_fit_casertano.pyx'], From 04baff1b4033527f7b302d5d4c47e2e3df98e3a7 Mon Sep 17 00:00:00 2001 From: PaulHuwe Date: Wed, 10 Jul 2024 12:43:10 -0400 Subject: [PATCH 7/9] Idea --- romanisim/tests/test_l3.py | 28 +++++++++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) diff --git a/romanisim/tests/test_l3.py b/romanisim/tests/test_l3.py index 34a7838c..51c31e10 100644 --- a/romanisim/tests/test_l3.py +++ b/romanisim/tests/test_l3.py @@ -74,6 +74,7 @@ def test_inject_sources_into_mosaic(): l3_mos.var_poisson.value[0:100, 100:200] = 0.02**2 l3_mos.var_poisson.value[100:200, 0:100] = 0.05**2 l3_mos.var_poisson.value[100:200, 100:200] = 0.1**2 + # l3_mos.var_poisson /= unit_factor.value**2 # Create normalized psf source catalog (same source in each quadrant) sc_dict = {"ra": 4 * [0.0], "dec": 4 * [0.0], "type": 4 * ["PSF"], "n": 4 * [-1.0], @@ -112,6 +113,7 @@ def test_inject_sources_into_mosaic(): # Set new poisson variance l3_mos.var_poisson = (l3_mos.data.value / Ct_all) * l3_mos.var_poisson.unit + # l3_mos.var_poisson = (l3_mos.data.value / Ct_all) * (l3_mos.var_poisson.unit / unit_factor.value**2) # Ensure that every data pixel value has increased or # remained the same with the new sources injected @@ -120,10 +122,34 @@ def test_inject_sources_into_mosaic(): # Ensure that every pixel's poisson variance has increased or # remained the same with the new sources injected # Numpy isclose is needed to determine equality, due to float precision issues - close_mask = np.isclose(l3_mos.var_poisson.value, l3_mos_orig.var_poisson.value, rtol=1e-06) + close_mask = np.isclose(l3_mos.var_poisson.value, l3_mos_orig.var_poisson.value, rtol=1e-20) + + import plotly.express as px + fig1 = px.imshow(l3_mos_orig.data.value, title='Mosaic Data', labels={'color':'counts / second', 'x':'y axis', 'y':'x axis'}) + fig1.show() + + fig2 = px.imshow(l3_mos.data.value, title='Sources Data', labels={'color':'counts / second', 'x':'y axis', 'y':'x axis'}) + fig2.show() + + fig3 = px.imshow((l3_mos.data.value - l3_mos_orig.data.value), title='Diff Data', labels={'color':'counts / second', 'x':'y axis', 'y':'x axis'}) + fig3.show() + + fig4 = px.imshow(close_mask, title='Mask', labels={'color': 'True / False', 'x':'y axis', 'y':'x axis'}) + fig4.show() + + fig5 = px.imshow(l3_mos_orig.var_poisson.value, title='Mosaic Poisson', labels={'color':'counts / second', 'x':'y axis', 'y':'x axis'}) + fig5.show() + + fig6 = px.imshow(l3_mos.var_poisson.value, title='Sources Poisson', labels={'color':'counts / second', 'x':'y axis', 'y':'x axis'}) + fig6.show() + + fig7 = px.imshow((l3_mos.var_poisson.value - l3_mos_orig.var_poisson.value), title='Diff Poisson', labels={'color':'counts / second', 'x':'y axis', 'y':'x axis'}) + fig7.show() + assert False in close_mask assert np.all(l3_mos.var_poisson.value[~close_mask] > l3_mos_orig.var_poisson.value[~close_mask]) + assert np.all(l3_mos.var_poisson.value > l3_mos_orig.var_poisson.value) # Create log entry and artifacts log.info('DMS232 successfully injected sources into a mosaic at points (50,50), (50,150), (150,50), (150,150).') From 657c83b7efb1f412786a7c2ac4d8b1767d2d8742 Mon Sep 17 00:00:00 2001 From: PaulHuwe Date: Wed, 24 Jul 2024 16:40:03 -0400 Subject: [PATCH 8/9] Integrated changes from L3 integration as well as addressed PR comments. --- romanisim/image.py | 25 ++++++++-------------- romanisim/l3.py | 30 +++++++++++++++----------- romanisim/parameters.py | 30 +++++++++++++++++++------- romanisim/tests/test_l3.py | 44 +++++++++++--------------------------- romanisim/wcs.py | 2 +- 5 files changed, 61 insertions(+), 70 deletions(-) diff --git a/romanisim/image.py b/romanisim/image.py index 9bc062e1..e350a3b6 100644 --- a/romanisim/image.py +++ b/romanisim/image.py @@ -214,9 +214,9 @@ def trim_objlist(objlist, image): def add_objects_to_image(image, objlist, xpos, ypos, psf, - flux_to_counts_factor, exptimes=None, + flux_to_counts_factor, convtimes=None, bandpass=None, filter_name=None, - wcs=None, rng=None, seed=None): + rng=None, seed=None): """Add sources to an image. Note: this includes Poisson noise when photon shooting is used @@ -236,16 +236,14 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf, physical fluxes in objlist (whether in profile SEDs or flux arrays) should be multiplied by this factor to convert to total counts in the image - exptimes: array_like - Exposure times to scale back to rate units + convtimes: array_like + Exposure times with unit scaling to convert to output rate units bandpass : galsim.Bandpass bandpass in which image is being rendered. This is used only in cases where chromatic profiles & PSFs are being used. filter_name : str filter to use to select appropriate flux from objlist. This is only used when achromatic PSFs and sources are being rendered. - wcs : galsim.GSFitsWCS - WCS corresponding to image rng : galsim.BaseDeviate random number generator to use seed : int @@ -281,10 +279,7 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf, for i, obj in enumerate(objlist): t0 = time.time() image_pos = galsim.PositionD(xpos[i], ypos[i]) - if wcs is None: - pwcs = image.wcs.local(image_pos) - else: - pwcs = wcs.local(image_pos) + pwcs = image.wcs.local(image_pos) profile = obj.profile if not chromatic: if obj.flux is None: @@ -295,10 +290,8 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf, psf0 = psf.at_position(xpos[i], ypos[i]) else: psf0 = psf - if isinstance(flux_to_counts_factor, list): - final = galsim.Convolve(profile * flux_to_counts_factor[i], psf0) - else: - final = galsim.Convolve(profile * flux_to_counts_factor, psf0) + factor = flux_to_counts_factor[i] if isinstance(flux_to_counts_factor, list) else flux_to_counts_factor + final = galsim.Convolve(profile * factor, psf0) if chromatic: stamp = final.drawImage( bandpass, center=image_pos, wcs=pwcs, @@ -310,8 +303,8 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf, except galsim.GalSimFFTSizeError: log.warning(f'Skipping source {i} due to too ' f'large FFT needed for desired accuracy.') - if exptimes is not None: - stamp /= exptimes[i] + if convtimes is not None: + stamp /= convtimes[i] bounds = stamp.bounds & image.bounds if bounds.area() > 0: diff --git a/romanisim/l3.py b/romanisim/l3.py index f6f0bbd8..9b988a3b 100644 --- a/romanisim/l3.py +++ b/romanisim/l3.py @@ -21,18 +21,17 @@ import romanisim.persistence from romanisim import log import roman_datamodels.maker_utils as maker_utils - -# Define centermost SCA for PSFs -CENTER_SCA = 2 +import roman_datamodels.datamodels as rdm +from roman_datamodels.stnode import WfiMosaic def add_objects_to_l3(l3_mos, source_cat, exptimes, xpos=None, ypos=None, coords=None, unit_factor=1.0, - coords_unit='rad', wcs=None, psf=None, rng=None, seed=None): + filter_name=None, coords_unit='rad', wcs=None, psf=None, rng=None, seed=None): """Add objects to a Level 3 mosaic Parameters ---------- - l3_mos : MosaicModel + l3_mos : MosaicModel or galsim.Image Mosaic of images source_cat : list List of catalog objects to add to l3_mos @@ -60,9 +59,9 @@ def add_objects_to_l3(l3_mos, source_cat, exptimes, xpos=None, ypos=None, coords None l3_mos is updated in place """ - # Obtain optical element - filter_name = l3_mos.meta.basic.optical_element + if filter_name is None: + filter_name = l3_mos.meta.basic.optical_element # Generate WCS (if needed) if wcs is None: @@ -70,10 +69,13 @@ def add_objects_to_l3(l3_mos, source_cat, exptimes, xpos=None, ypos=None, coords # Create PSF (if needed) if psf is None: - psf = romanisim.psf.make_psf(filter_name=filter_name, sca=CENTER_SCA, chromatic=False, webbpsf=True) + psf = romanisim.psf.make_psf(filter_name=filter_name, sca=parameters.default_sca, chromatic=False, webbpsf=True) # Create Image canvas to add objects to - sourcecountsall = galsim.ImageF(l3_mos.data.value, wcs=wcs, xmin=0, ymin=0) + if isinstance(l3_mos, (rdm.MosaicModel, WfiMosaic)): + sourcecountsall = galsim.ImageF(l3_mos.data.value, wcs=wcs, xmin=0, ymin=0) + else: + sourcecountsall = l3_mos # Create position arrays (if needed) if any(pos is None for pos in [xpos, ypos]): @@ -87,12 +89,14 @@ def add_objects_to_l3(l3_mos, source_cat, exptimes, xpos=None, ypos=None, coords # Add sources to the original mosaic data array romanisim.image.add_objects_to_image(sourcecountsall, source_cat, xpos=xpos, ypos=ypos, - psf=psf, flux_to_counts_factor=[xpt * unit_factor for xpt in exptimes], - exptimes=exptimes, bandpass=[filter_name], filter_name=filter_name, - wcs=wcs, rng=rng, seed=seed) + psf=psf, flux_to_counts_factor=exptimes, + convtimes=[xpt / unit_factor for xpt in exptimes], + bandpass=[filter_name], filter_name=filter_name, + rng=rng, seed=seed) # Save array with added sources - l3_mos.data = sourcecountsall.array * l3_mos.data.unit + if isinstance(l3_mos, (rdm.MosaicModel, WfiMosaic)): + l3_mos.data = sourcecountsall.array * l3_mos.data.unit # l3_mos is updated in place, so no return return None diff --git a/romanisim/parameters.py b/romanisim/parameters.py index 0003eb74..016a3932 100644 --- a/romanisim/parameters.py +++ b/romanisim/parameters.py @@ -84,14 +84,25 @@ "linearity": None, "readnoise": 5.0 * u.DN, "saturation": 55000 * u.DN, - "photom": {"F062": 0.28798833 * u.MJy / u.sr, - "F087": 0.3911463 * u.MJy / u.sr, - "F106": 0.37039083 * u.MJy / u.sr, - "F129": 0.36495176 * u.MJy / u.sr, - "F146": 0.11750617 * u.MJy / u.sr, - "F158": 0.35352992 * u.MJy / u.sr, - "F184": 0.53612772 * u.MJy / u.sr, - "F213": 0.55763922 * u.MJy / u.sr, + # Taken from roman_wfi_photom_0046.asdf + "photom": {"photmjsr": {"F062": 0.28798833 * u.MJy / u.sr, + "F087": 0.3911463 * u.MJy / u.sr, + "F106": 0.37039083 * u.MJy / u.sr, + "F129": 0.36495176 * u.MJy / u.sr, + "F146": 0.11750617 * u.MJy / u.sr, + "F158": 0.35352992 * u.MJy / u.sr, + "F184": 0.53612772 * u.MJy / u.sr, + "F213": 0.55763922 * u.MJy / u.sr, + }, + "pixelareasr": {"F062": 2.640600009241359e-13 * u.sr, + "F087": 2.640600009241359e-13 * u.sr, + "F106": 2.640600009241359e-13 * u.sr, + "F129": 2.640600009241359e-13 * u.sr, + "F146": 2.640600009241359e-13 * u.sr, + "F158": 2.640600009241359e-13 * u.sr, + "F184": 2.640600009241359e-13 * u.sr, + "F213": 2.640600009241359e-13 * u.sr, + } } } @@ -149,6 +160,9 @@ "pixel_depth": 5, } +# Centermost PSF to use for mosaic creation +default_sca = 2 + # Persistence defaults # delete persistence records fainter than 0.01 electron / s # e.g., MA table 1 has ~144 s, so this is ~1 electron over the whole exposure. diff --git a/romanisim/tests/test_l3.py b/romanisim/tests/test_l3.py index 51c31e10..c1a36139 100644 --- a/romanisim/tests/test_l3.py +++ b/romanisim/tests/test_l3.py @@ -34,7 +34,8 @@ def test_inject_sources_into_mosaic(): galsim.roman.n_pix = 200 rng_seed = 42 metadata = copy.deepcopy(parameters.default_mosaic_parameters_dictionary) - metadata['basic']['optical_element'] = 'F158' + filter_name = 'F158' + metadata['basic']['optical_element'] = filter_name # Create WCS twcs = wcs.get_mosaic_wcs(metadata, shape=(galsim.roman.n_pix, galsim.roman.n_pix)) @@ -53,6 +54,7 @@ def test_inject_sources_into_mosaic(): # Create level 3 mosaic model l3_mos = maker_utils.mk_level3_mosaic(shape=(galsim.roman.n_pix, galsim.roman.n_pix)) + l3_mos['meta']['wcs'] = twcs # Update metadata in the l3 model for key in metadata.keys(): @@ -60,7 +62,9 @@ def test_inject_sources_into_mosaic(): l3_mos.meta[key].update(metadata[key]) # Obtain unit conversion factor - unit_factor = parameters.reference_data['photom'][l3_mos.meta.basic.optical_element] + # Need to convert from counts / pixel to MJy / sr + unit_factor = ((3631 * u.Jy) / (romanisim.bandpass.get_abflux(filter_name) + * parameters.reference_data['photom']["pixelareasr"][filter_name])).to(u.MJy / u.sr) # Populate the mosaic data array with gaussian noise from generators g1.generate(l3_mos.data.value[0:100, 0:100]) @@ -74,11 +78,11 @@ def test_inject_sources_into_mosaic(): l3_mos.var_poisson.value[0:100, 100:200] = 0.02**2 l3_mos.var_poisson.value[100:200, 0:100] = 0.05**2 l3_mos.var_poisson.value[100:200, 100:200] = 0.1**2 - # l3_mos.var_poisson /= unit_factor.value**2 + l3_mos.var_poisson *= unit_factor.value**2 # Create normalized psf source catalog (same source in each quadrant) sc_dict = {"ra": 4 * [0.0], "dec": 4 * [0.0], "type": 4 * ["PSF"], "n": 4 * [-1.0], - "half_light_radius": 4 * [0.0], "pa": 4 * [0.0], "ba": 4 * [1.0], "F158": 4 * [1.0]} + "half_light_radius": 4 * [0.0], "pa": 4 * [0.0], "ba": 4 * [1.0], filter_name: 4 * [1.0]} sc_table = table.Table(sc_dict) # Set locations @@ -97,7 +101,7 @@ def test_inject_sources_into_mosaic(): sc_table["ra"][idx], sc_table["dec"][idx] = (twcs._radec(x, y) * u.rad).to(u.deg).value - source_cat = catalog.table_to_catalog(sc_table, ["F158"]) + source_cat = catalog.table_to_catalog(sc_table, [filter_name]) # Copy original Mosaic before adding sources as sources are added in place l3_mos_orig = l3_mos.copy() @@ -109,11 +113,11 @@ def test_inject_sources_into_mosaic(): # Create overall scaling factor map Ct_all = np.divide(l3_mos_orig.data.value, l3_mos_orig.var_poisson.value, - out=np.ones(l3_mos_orig.data.shape), where=l3_mos_orig.var_poisson.value != 0) + out=np.ones(l3_mos_orig.data.shape), + where=l3_mos_orig.var_poisson.value != 0) # Set new poisson variance l3_mos.var_poisson = (l3_mos.data.value / Ct_all) * l3_mos.var_poisson.unit - # l3_mos.var_poisson = (l3_mos.data.value / Ct_all) * (l3_mos.var_poisson.unit / unit_factor.value**2) # Ensure that every data pixel value has increased or # remained the same with the new sources injected @@ -122,34 +126,10 @@ def test_inject_sources_into_mosaic(): # Ensure that every pixel's poisson variance has increased or # remained the same with the new sources injected # Numpy isclose is needed to determine equality, due to float precision issues - close_mask = np.isclose(l3_mos.var_poisson.value, l3_mos_orig.var_poisson.value, rtol=1e-20) - - import plotly.express as px - fig1 = px.imshow(l3_mos_orig.data.value, title='Mosaic Data', labels={'color':'counts / second', 'x':'y axis', 'y':'x axis'}) - fig1.show() - - fig2 = px.imshow(l3_mos.data.value, title='Sources Data', labels={'color':'counts / second', 'x':'y axis', 'y':'x axis'}) - fig2.show() - - fig3 = px.imshow((l3_mos.data.value - l3_mos_orig.data.value), title='Diff Data', labels={'color':'counts / second', 'x':'y axis', 'y':'x axis'}) - fig3.show() - - fig4 = px.imshow(close_mask, title='Mask', labels={'color': 'True / False', 'x':'y axis', 'y':'x axis'}) - fig4.show() - - fig5 = px.imshow(l3_mos_orig.var_poisson.value, title='Mosaic Poisson', labels={'color':'counts / second', 'x':'y axis', 'y':'x axis'}) - fig5.show() - - fig6 = px.imshow(l3_mos.var_poisson.value, title='Sources Poisson', labels={'color':'counts / second', 'x':'y axis', 'y':'x axis'}) - fig6.show() - - fig7 = px.imshow((l3_mos.var_poisson.value - l3_mos_orig.var_poisson.value), title='Diff Poisson', labels={'color':'counts / second', 'x':'y axis', 'y':'x axis'}) - fig7.show() - + close_mask = np.isclose(l3_mos.var_poisson.value, l3_mos_orig.var_poisson.value, rtol=1e-06) assert False in close_mask assert np.all(l3_mos.var_poisson.value[~close_mask] > l3_mos_orig.var_poisson.value[~close_mask]) - assert np.all(l3_mos.var_poisson.value > l3_mos_orig.var_poisson.value) # Create log entry and artifacts log.info('DMS232 successfully injected sources into a mosaic at points (50,50), (50,150), (150,50), (150,150).') diff --git a/romanisim/wcs.py b/romanisim/wcs.py index 59c31380..7a91293d 100644 --- a/romanisim/wcs.py +++ b/romanisim/wcs.py @@ -441,7 +441,7 @@ def get_mosaic_wcs(mosaic, shape=None): mosaic_node['meta'][key].update(mosaic[key]) else: mosaic_node['meta'][key] = mosaic[key] - mosaic_node = roman_datamodels.datamodels.MosaicModel(mosaic_node) + # mosaic_node = roman_datamodels.datamodels.MosaicModel(mosaic_node) else: mosaic_node = mosaic From d4af8fd932768c45a9fd5e0d6938bc3c3a3240c0 Mon Sep 17 00:00:00 2001 From: PaulHuwe Date: Tue, 6 Aug 2024 15:47:25 -0400 Subject: [PATCH 9/9] Removed extraneous comment. --- romanisim/wcs.py | 1 - 1 file changed, 1 deletion(-) diff --git a/romanisim/wcs.py b/romanisim/wcs.py index 7a91293d..c0f7c289 100644 --- a/romanisim/wcs.py +++ b/romanisim/wcs.py @@ -460,6 +460,5 @@ def get_mosaic_wcs(mosaic, shape=None): world_origin=galsim.PositionD(0, 0)) wcs = galsim.TanWCS(affine, util.celestialcoord(world_pos)) - # wcs = GWCS(wcs) return wcs