diff --git a/docs/index.rst b/docs/index.rst index 85bc4944..683294f5 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -19,6 +19,7 @@ Contents romanisim/image romanisim/l1 romanisim/l2 + romanisim/l3 romanisim/refs romanisim/catalog diff --git a/docs/romanisim/image.rst b/docs/romanisim/image.rst index 09a399dc..4df992fc 100644 --- a/docs/romanisim/image.rst +++ b/docs/romanisim/image.rst @@ -14,4 +14,19 @@ The parsed metadata is used to make a ``counts`` image that is an idealized imag These idealized count images are then used to either make a level 2 image or a level 1 image, which are intended to include the effects of these complications. The construction of L1 images is described :doc:`here `, and the construction of L2 images is described :doc:`here `. +L2 source injection +------------------- + +We support injecting sources into L2 files. L2 source injection follows many of the same steps as L1 creation, with some short cuts. The simulation creates an idealized image of the total number of counts deposited into each pixel, using the provided catalog, exposure time, filter, WCS, and PSF from the provided image. A virtual ramp is generated for the existing L2 file by evenly apportioning the L2 flux among the resultants. Additional counts from the injected source are added to each resultant in the virtual ramp using the same code as for the L1 simulations. The resulting virtual ramp is fit to get the new flux per pixel, and we replace the values in the original L2 model with the new slope measurements. + +This makes some shortcuts: + +* The L2 files don't include information about which pixels contain CR hits. Sources injected into pixels with CR hits get re-fit as if there were no CR. +* Non-linearity is not included; the ramp is refit without any non-linearity. +* The ramp is refit without persistence, but any persistence which was included in the initial ramp slope is still included. + +The function :meth:`romanisim.image.inject_sources_into_l2` is the intended entry point for L2 source injection; see its documentation for more details about how to add sources to an L2 image. + + + .. automodapi:: romanisim.image diff --git a/docs/romanisim/l2.rst b/docs/romanisim/l2.rst index e3c29223..e210699f 100644 --- a/docs/romanisim/l2.rst +++ b/docs/romanisim/l2.rst @@ -16,4 +16,16 @@ This is a fairly faithful representation of how level two image construction wor * We have a simplistic saturation treatment, just clipping resultants that exceed the saturation level to the saturation level and throwing a flag. + +L2 source injection +=================== +We support L2 source injection through :meth:`romanisim.image.inject_sources_into_l2`. This routine is intended for cases when existing L2 images are available, and a small number of sources need to be added on top of those images. The L2 source injection is intended to be a fairly faithful simulation of the process. It includes the following steps: +* The model WCS is used to generate the pixel locations for sources. +* The number of additional electrons in each pixel is simulated. +* A virtual ramp is created that evenly apportions the original L2 flux along the ramp. +* The new counts are apportioned randomly along the ramp following a Poisson distribution in the same way as for the usual L1 simulations. +* The new ramp is then fit using the ramp fitting code, and the new fit replaces the old fit in the L2 file. + +Note that this procedure does not include any new cosmic rays, or any special treatment of pixels that had cosmic rays in them and now have additional source flux. It is hard to do this right since we do not know which pixel the cosmic ray landed in on the basis of the L2 file. If you need more accuracy than this you probably need to directly inject sources into the L1 file. + .. automodapi:: romanisim.ramp diff --git a/docs/romanisim/l3.rst b/docs/romanisim/l3.rst new file mode 100644 index 00000000..bcd28fdc --- /dev/null +++ b/docs/romanisim/l3.rst @@ -0,0 +1,59 @@ +L3 images +========= + +We also simulate L3 images, also known as coadded images or mosaics; these terms are used equivalently here. In contrast to L1 and L2 images, L3 images do not do a full simulation of the individual reads entering each resultant and the following ramp fitting. Accordingly, they also do not include most systematic effects---for example, cosmic rays and non-linearity are not simulated. Instead, the L3 simulations conceptually are more closely related to the initial idealized :doc:`count images ` than they are to L1 or L2 images. They differ from these idealized images in that we attempt to change the PSF in order to account for non-native pixel scales that are possible in coadded images. + +The simulation of L3 images is simplified relative to true L3 images. First, as above, we do not do the full L2 simulations. Second, we do not have a notion of a series of images that are being coadded---we are instead just saying that one could produce a mosaic that has a particular pixel scale, effective read noise, and exposure time. We do not attempt to simulate the drizzle process used for typical L3 images in Roman, and accordingly do not include any drizzle-induced covariances among the pixels. Similarly, because we don't really have a notion of a series of input L2 images, the L3 PSF is not a combination of the PSF of the input images, but is instead taken as the PSF evaluated at the center of the WFI02 detector, and then additionally convolved depending on the output pixel grid scale. + +The L3 PSF aims to include the convolution with the native scale that Roman necessarily includes, but to also support oversampled images where the L3 pixel scale can be significantly finer than the native pixel scale. In these cases we include an additional square convolution to "make up the difference" between the L3 pixel scale and the L2 pixel scale; see the :meth:`romanisim.l3.l3_psf` method for more details. + +The ``romanisim-make-l3`` command-line interface can be used to make L3 files. Invocations look like:: + + romanisim-make-l3 --bandpass F158 --radec 270 66 --npix 2000 --pixscalefrac 0.9 --exptime 1000 l3.asdf gaia-270-66-2027-06-01.ecsv + +which would make an F158 image with 2000x2000 pixels, slightly oversampled, corresponding to 1000 seconds, using the given catalog and output filename. + +There are a number of possible input arguments:: + + romanisim-make-l3 -h + usage: romanisim-make-l3 [-h] [--bandpass BANDPASS] [--config CONFIG] [--date DATE] [--nobj NOBJ] [--radec RADEC RADEC] [--npix NPIX] + [--pixscalefrac PIXSCALEFRAC] [--exptime EXPTIME] [--effreadnoise EFFREADNOISE] [--nexposures NEXPOSURES] + [--rng_seed RNG_SEED] + filename catalog + + Make an L3 image. + + positional arguments: + filename output image (asdf) + catalog input catalog (ecsv) + + options: + -h, --help show this help message and exit + --bandpass BANDPASS bandpass to simulate (default: F087) + --config CONFIG input parameter override file (yaml) (default: None) + --date DATE UTC Date and Time of observation to simulate in ISOT format. (default: None) + --radec RADEC RADEC ra and dec (deg) (default: None) + --npix NPIX number of pixels across image (default: 4000) + --pixscalefrac PIXSCALEFRAC + pixel scale as fraction of original (default: 1.0) + --exptime EXPTIME total exposure time on field; roughly time per exposure times number of exposures (default: 100.0) + --effreadnoise EFFREADNOISE + effective readnoise per pixel in MJy/sr. If None, a pessimistic estimate is computed. (default: None) + --nexposures NEXPOSURES + number of exposures on field. Used only to compute the effective read noise. (default: 1) + --rng_seed RNG_SEED + + EXAMPLE: romanisim-make-l3 output_image.asdf catalog.ecsv + +The arguments control things like the filter used, the date of observation (used to estimate the sky background), the right ascension and declination of the center of the mosaic, the number of pixels in the mosaic, and the degree of oversampling of the output mosaic. Parameters also control the effective depth of the mosaic, through the exposure time, effective read noise, and number of exposures. + +L3 source injection +------------------- + +We support injecting sources into L3 files. L3 source injection is very similar to L3 image creation, except we do not need to invent effective read noises or exposure times or WCSes. Instead we only change the Poisson noise, since we are only adding additional new sources. + +Conceptually we just add more sources to the image in the same way as in L3 image creation. The primary challenge is picking the correct exposure time so that the Poisson noise in the simulation is correct; we (optionally) simulate individual photons and so the conversion between calibrated units and electrons is still relevant in the L3 files. We do this by using the relationship between the calibrated flux and corresponding Poisson variance to figure out the effective calibrated-flux-to-electron relationship, and then simulate the appropriate number of electrons and convert back. + +The function :meth:`romanisim.l3.inject_sources_into_l3` is the intended entry point for L3 source injection; see its documentation for more details about how to add sources to an L3 image. + +.. automodapi:: romanisim.l3 diff --git a/romanisim/bandpass.py b/romanisim/bandpass.py index d146b613..bbdb0ec6 100644 --- a/romanisim/bandpass.py +++ b/romanisim/bandpass.py @@ -12,6 +12,7 @@ from astropy.table import Table from astropy import units as u import functools +from romanisim import parameters # to go from calibrated fluxes in maggies to counts in the Roman bands # we need to multiply by a constant determined by the AB magnitude @@ -181,3 +182,26 @@ def compute_count_rate(flux, bandpass, filename=None, effarea=None, wavedist=Non # per second return zpflux + + +def etomjysr(bandpass): + """Compute factor converting e/s/pix to MJy/sr. + + Assumes a pixel scale of 0.11" (romanisim.parameters.pixel_scale) + + Parameters + ---------- + bandpass : str + the name of the bandpass + + Returns + ------- + float + the factor F such that MJy / sr = F * DN/s + """ + + abflux = get_abflux(bandpass) # e/s corresponding to 3631 Jy + srpix = ((parameters.pixel_scale * u.arcsec) ** 2).to(u.sr).value + mjysr = 1 / abflux * 3631 / 10 ** 6 / srpix # should be ~0.3 + return mjysr + diff --git a/romanisim/image.py b/romanisim/image.py index b9ded592..20be1e32 100644 --- a/romanisim/image.py +++ b/romanisim/image.py @@ -205,7 +205,7 @@ def trim_objlist(objlist, image): def add_objects_to_image(image, objlist, xpos, ypos, psf, - flux_to_counts_factor, convtimes=None, + flux_to_counts_factor, outputunit_to_electrons=None, bandpass=None, filter_name=None, add_noise=False, rng=None, seed=None): """Add sources to an image. @@ -227,8 +227,9 @@ 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 - convtimes: array_like - Exposure times with unit scaling to convert to output rate units + outputunit_to_electrons : array_like + One output image unit corresponds to this many electrons. If None, + leave as electrons. bandpass : galsim.Bandpass bandpass in which image is being rendered. This is used only in cases where chromatic profiles & PSFs are being used. @@ -287,17 +288,17 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf, stamp = final.drawImage( bandpass, center=image_pos, wcs=pwcs, method='phot', rng=rng) - if add_noise: - stamp.addNoise(galsim.PoissonNoise(rng)) else: try: stamp = final.drawImage(center=image_pos, wcs=pwcs) + if add_noise: + stamp.addNoise(galsim.PoissonNoise(rng)) except galsim.GalSimFFTSizeError: log.warning(f'Skipping source {i} due to too ' f'large FFT needed for desired accuracy.') - if convtimes is not None: - stamp /= convtimes[i] + if outputunit_to_electrons is not None: + stamp /= outputunit_to_electrons[i] bounds = stamp.bounds & image.bounds if bounds.area() > 0: @@ -599,7 +600,8 @@ def gather_reference_data(image_mod, usecrds=False): """ reffiles = {k: v for k, v in parameters.reference_data.items()} - reffiles.pop('photom') + if 'photom' in reffiles: + reffiles.pop('photom') out = dict(**reffiles) if usecrds: @@ -917,7 +919,7 @@ def inject_sources_into_l2(model, cat, x=None, y=None, psf=None, rng=None, This routine allows sources to be injected onto an existing L2 image. Source injection into an L2 image relies on knowing the objects' x and y locations, the PSF, and the image gain; if these are not provided, - reasonable defaults are generated. + reasonable defaults are generated from the input model. The simulation proceeds by (optionally) using the model WCS to generate the x & y locations, grabbing the gain from @@ -1036,7 +1038,7 @@ def inject_sources_into_l2(model, cat, x=None, y=None, psf=None, rng=None, newramp[:, m], read_pattern, gain=gain, flat=1, darkrate=0) - res = model.copy() + res = copy.deepcopy(model) res.data[m] = newimage res.var_rnoise[m] = readvar res.var_poisson[m] = poissonvar diff --git a/romanisim/l3.py b/romanisim/l3.py index c7708d53..a793af68 100644 --- a/romanisim/l3.py +++ b/romanisim/l3.py @@ -24,40 +24,37 @@ import astropy from galsim import roman from astropy import table +import astropy.coordinates -def add_objects_to_l3(l3_mos, source_cat, exptimes, xpos=None, ypos=None, coords=None, cps_conv=1.0, unit_factor=1.0, - filter_name=None, bandpass=None, coords_unit='rad', wcs=None, psf=None, rng=None, seed=None): +def add_objects_to_l3(l3_mos, source_cat, exptimes, xpos, ypos, psf, + etomjysr=1.0, maggytoes=1.0, + filter_name=None, bandpass=None, + rng=None, seed=None): """Add objects to a Level 3 mosaic Parameters ---------- l3_mos : MosaicModel or galsim.Image Mosaic of images - source_cat : list + source_cat : list[CatalogObject] 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 - cps_conv : float - Factor to convert data to cps - unit_factor: float - Factor to convert counts data to MJy / sr + psf : galsim.Profile + PSF to use + etomjysr: float + Factor to convert electrons to MJy / sr + maggytoes : float + Factor to convert maggies to e/s filter_name : str Filter to use to select appropriate flux from objlist. This is only used when achromatic PSFs and sources are being rendered. bandpass : galsim.Bandpass Bandpass in which mosaic is being rendered. This is used only in cases where chromatic profiles & PSFs are being used. - 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 @@ -65,59 +62,129 @@ def add_objects_to_l3(l3_mos, source_cat, exptimes, xpos=None, ypos=None, coords Returns ------- - None - l3_mos is updated in place + Information from romanisim.image.add_objects_to_image. Note + that l3_mos is updated in place. """ # Obtain optical element if filter_name is None: filter_name = l3_mos.meta.basic.optical_element - if bandpass is None: - bandpass = [filter_name] - - # Generate WCS (if needed) - if wcs is None: - wcs = romanisim.wcs.get_mosaic_wcs(l3_mos.meta, shape=l3_mos.data.shape) - - # Create PSF (if needed) - if psf is None: - psf = romanisim.psf.make_psf(filter_name=filter_name, sca=romanisim.parameters.default_sca, chromatic=False, webbpsf=True) # Create Image canvas to add objects to if isinstance(l3_mos, (rdm.MosaicModel, WfiMosaic)): - sourcecountsall = galsim.ImageF(l3_mos.data.value, wcs=wcs, xmin=0, ymin=0) + sourcecountsall = galsim.ImageF( + l3_mos.data.value, wcs=romanisim.wcs.GWCS(l3_mos.meta.wcs), + xmin=0, ymin=0) else: sourcecountsall = l3_mos - # 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' - # Generate x,y positions for sources - xpos, ypos = sourcecountsall.wcs.radecToxy(coords[:, 0], coords[:, 1], coords_unit) - # Add sources to the original mosaic data array - outinfo = romanisim.image.add_objects_to_image(sourcecountsall, source_cat, xpos=xpos, ypos=ypos, - psf=psf, flux_to_counts_factor=[xpt * cps_conv for xpt in exptimes], - convtimes=[xpt / unit_factor for xpt in exptimes], - bandpass=bandpass, filter_name=filter_name, - rng=rng, seed=seed) + outinfo = romanisim.image.add_objects_to_image( + sourcecountsall, source_cat, xpos, ypos, + psf, flux_to_counts_factor=[xpt * maggytoes for xpt in exptimes], + outputunit_to_electrons=[xpt / etomjysr for xpt in exptimes], + bandpass=bandpass, filter_name=filter_name, + rng=rng, seed=seed, add_noise=True) # Save array with added sources 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 outinfo -def generate_mosaic_geometry(): - """Generate a geometry map (context) for a mosaic dither pattern / set of pointings and roll angles - TBD +def inject_sources_into_l3(model, cat, x=None, y=None, psf=None, rng=None, + webbpsf=True, exptimes=None, seed=None): + """Inject sources into an L3 image. + + This routine allows sources to be injected onto an existing L3 image. + Source injection into an L3 image relies on knowing the objects' + x and y locations, the PSF, and the exposure time; if these are not + provided, reasonable defaults are generated. + + The simulation proceeds by (optionally) using the model WCS to generate the + x & y locations. It also optionally computes an appropriate pixel-convolved + PSF for the image. It optionally uses the model Poisson variances and + model fluxes to estimate the per-source exposure times appropriate + for injecting the sources in the catalog onto the image. Finally, + it updates the var_poisson of the input image to account for the additional + variance from the added sources. + + Parameters + ---------- + model: roman_datamodels.datamodels.WfiMosaic + model into which to inject sources + cat: astropy.table.Table + catalog of sources to inject into image + x: list[float] or None + x coordinates of catalog locations in image + y: list[float] or None + y coordinates of catalog locations in image + exptimes : list[float] or None + exposure times of each source. Computed from model.var_poisson + and model.flux at each source location if not specified. + psf: galsim.gsobject.GSObject + PSF to use + rng: galsim.BaseDeviate + galsim random number generator to use + seed : int + Seed to use for rng + webbpsf: bool + if True, use WebbPSF to model the PSF + + Returns + ------- + outinfo : np.ndarray with information about added sources """ - pass + if seed is None: + seed = 125 + if rng is None: + rng = galsim.UniformDeviate(seed) + + if x is None or y is None: + x, y = model.meta.wcs.numerical_inverse(cat['ra'], cat['dec'], + with_bounding_box=False) + + filter_name = model.meta.basic.optical_element + cat = romanisim.catalog.table_to_catalog(cat, [filter_name]) + + wcs = romanisim.wcs.GWCS(model.meta.wcs) + pixscalefrac = get_pixscalefrac(wcs, model.data.shape) + if psf is None: + if (pixscalefrac > 1) or (pixscalefrac < 0): + raise ValueError('weird pixscale!') + psf = l3_psf(filter_name, pixscalefrac, webbpsf=True, chromatic=False) + + maggytoes = romanisim.bandpass.get_abflux(filter_name) + etomjysr = romanisim.bandpass.etomjysr(filter_name) / pixscalefrac ** 2 + + Ct = [] + for idx, (x0, y0) in enumerate(zip(x, y)): + # Set scaling factor for injected sources + # Flux / sigma_p^2 + xidx, yidx = int(np.round(x0)), int(np.round(y0)) + if model.var_poisson[yidx, xidx].value != 0: + Ct.append(math.fabs( + model.data[yidx, xidx].value / + model.var_poisson[yidx, xidx].value)) + else: + Ct.append(1.0) + Ct = np.array(Ct) + # etomjysr = 1/C; C converts fluxes to electrons + exptimes = Ct * etomjysr + + Ct_all = (model.data.value / + (model.var_poisson.value + (model.var_poisson.value == 0))) + + # compute the total number of counts we got from the source + res = add_objects_to_l3( + model, cat, exptimes, x, y, psf, etomjysr=etomjysr, + maggytoes=maggytoes, filter_name=filter_name, bandpass=None, + rng=rng) + + model.var_poisson = (model.data.value / Ct_all) * model.var_poisson.unit + + return res + def generate_exptime_array(cat, meta): @@ -127,7 +194,7 @@ def generate_exptime_array(cat, meta): """ # Get wcs for this metadata - twcs = romanisim.wcs.get_mosaic_wcs(meta) + twcs = romanisim.wcs.GWCS(romanisim.wcs.get_mosaic_wcs(meta)) # Obtain sky positions for objects coords = np.array([[o.sky_pos.ra.rad, o.sky_pos.dec.rad] @@ -155,10 +222,63 @@ def generate_exptime_array(cat, meta): return context -def simulate(shape, wcs, efftimes, filter, catalog, metadata={}, - effreadnoise=None, sky=None, psf=None, coords_unit='rad', - cps_conv=None, unit_factor=None, - bandpass=None, seed=None, rng=None, +def l3_psf(bandpass, scale=0, chromatic=False, **kw): + """Construct a PSF for an L3 image. + + The challenge here is that the L3 PSF has gone through drizzling which + broadens the PSF relative to the native PSF. We treat this here by doing + the following: + + * we start with the native PSF for the appropriate bandpass for SCA 2 + (no pixel convolution) + * we convolve that with a box of size sqrt(1-scale**2) pixel + * we return the result + + So for a very small scale you still get a convolution by the full + native pixel scaling (since the image was ultimately observed with the + Roman detectors and needs to be seen through that grid), but you get + oversampled output. There will + always be an additional convolution with the output pixel scale. In + the limit that the output pixel scale is 1, this does nothing, but + provides undersampled output. + + Extra arguments are passed to romanisim.psf.make_psf. + + Parameters + ---------- + scale : float + The output mosaic pixel scale. Must be between 0 and 1. + bandpass : str + The filter to use + chromatic : bool + if True, generate a chromatic PSF rather than an achromatic PSF. + This is intended for use with galsim chromatic sources. + + Returns + ------- + galsim.Profile object for use as the PSF + """ + + if scale < 0 or scale > 1: + raise ValueError('scale must be between 0 and 1') + convscale = romanisim.parameters.pixel_scale * np.sqrt( + 1 - scale**2) + if scale == 1: + extra_convolution = None + else: + extra_convolution = galsim.Pixel( + convscale * romanisim.parameters.pixel_scale) + psf = romanisim.psf.make_psf(filter_name=bandpass, + sca=romanisim.parameters.default_sca, + extra_convolution=extra_convolution, + chromatic=chromatic, **kw) + return psf + + +def simulate(shape, wcs, efftimes, filter_name, catalog, nexposures=1, + metadata={}, + effreadnoise=None, sky=None, psf=None, + bandpass=None, seed=None, rng=None, webbpsf=True, **kwargs): """Simulate a sequence of observations on a field in different bandpasses. @@ -166,35 +286,34 @@ def simulate(shape, wcs, efftimes, filter, catalog, metadata={}, ---------- shape : tuple Array shape of mosaic to simulate. - wcs : galsim.GSFitsWCS - WCS corresponding to image + wcs : gwcs.wcs.WCS + WCS corresponding to image. Will only work well with square pixels. efftimes : np.ndarray or float - Effective exposure time of reach pixel in mosaic. + Time Roman spent observing each part of the sky. If an array, shape must match shape parameter. - filter : str + filter_name : str Filter to use to select appropriate flux from objlist. This is only used when achromatic PSFs and sources are being rendered. catalog : list[CatalogObject] or Table List of catalog objects to add to l3_mos + nexposures : int + Number of exposures on the field. Used to compute effreadnoise + estimate if not explicitly set, otherwise not used. metadata : dict Metadata structure for Roman asdf file. effreadnoise : float - Effective read noise for mosaic. + Effective read noise for mosaic (MJy / sr) sky : float or array_like - Image or constant with the counts / pix / sec from sky. If None, then + Image or constant with sky and other backgrounds (MJy / sr). If None, then sky will be generated from galsim's getSkyLevel for Roman for the date provided in metadata[basic][time_mean_mjd]. - psf : galsim.Profile + psf : galsim.Profile or None PSF for image - coords_unit : string - units of coords - cps_conv : float - Factor to convert data to cps - unit_factor: float - Factor to convert counts data to MJy / sr bandpass : galsim.Bandpass Bandpass in which mosaic is being rendered. This is used only in cases where chromatic profiles & PSFs are being used. + webbpsf : bool + Use webbpsf to compute PSF rng : galsim.BaseDeviate random number generator to use seed : int @@ -213,26 +332,34 @@ def simulate(shape, wcs, efftimes, filter, catalog, metadata={}, # Create metadata object meta = maker_utils.mk_mosaic_meta() for key in romanisim.parameters.default_mosaic_parameters_dictionary.keys(): - meta[key].update(romanisim.parameters.default_mosaic_parameters_dictionary[key]) + meta[key].update( + romanisim.parameters.default_mosaic_parameters_dictionary[key]) meta['wcs'] = wcs - meta['basic']['optical_element'] = filter + meta['basic']['optical_element'] = filter_name for key in metadata.keys(): - if key not in meta: + if key not in meta or not isinstance(meta[key], dict): meta[key] = metadata[key] else: meta[key].update(metadata[key]) - # May need an equivalent of this this for L3? - # util.add_more_metadata(meta) - log.info('Simulating filter {0}...'.format(filter)) + add_more_metadata(meta, efftimes, filter_name, wcs, shape, nexposures) + + log.info('Simulating filter {0}...'.format(filter_name)) # Get filter and bandpass - galsim_filter_name = romanisim.bandpass.roman2galsim_bandpass[filter] + galsim_filter_name = romanisim.bandpass.roman2galsim_bandpass[filter_name] if bandpass is None: bandpass = roman.getBandpasses(AB_zeropoint=True)[galsim_filter_name] - # Create initial galsim image (x & y are flipped) - image = galsim.ImageF(shape[1], shape[0], wcs=wcs, xmin=0, ymin=0) + # Create initial galsim image; galsim wants x/y instead of normal shape + image = galsim.ImageF(shape[1], shape[0], wcs=romanisim.wcs.GWCS(wcs), + xmin=0, ymin=0) + + pixscalefrac = get_pixscalefrac(image.wcs, shape) + etomjysr = romanisim.bandpass.etomjysr(filter_name) / pixscalefrac ** 2 + # this should really be per-pixel to deal with small distortions, + # but these are 0.01% 1 degree away in a tangent plane projection, + # and we ignore them. # Create sky for this mosaic, if not provided (in cps) if sky is None: @@ -240,103 +367,133 @@ def simulate(shape, wcs, efftimes, filter, catalog, metadata={}, if not isinstance(date, astropy.time.Time): date = astropy.time.Time(date, format='mjd') - # Examine this in more detail to ensure it is correct - # Create base sky level - mos_cent_pos = wcs.toWorld(image.true_center) + mos_cent_pos = image.wcs.toWorld(image.true_center) sky_level = roman.getSkyLevel(bandpass, world_pos=mos_cent_pos, exptime=1) sky_level *= (1.0 + roman.stray_light_fraction) - sky_mosaic = image * 0 - wcs.makeSkyImage(sky_mosaic, sky_level) - sky_mosaic += roman.thermal_backgrounds[galsim_filter_name] - sky = sky_mosaic - - # Obtain unit conversion factors - # Need to convert from counts / pixel to MJy / sr - # Flux to counts - if cps_conv is None: - cps_conv = romanisim.bandpass.get_abflux(filter) - # Unit factor - if unit_factor is None: - unit_factor = ((3631 * u.Jy) / (romanisim.bandpass.get_abflux(filter) * 10e6 - * romanisim.parameters.reference_data['photom']["pixelareasr"][filter])).to(u.MJy / u.sr) + sky = image * 0 + image.wcs.makeSkyImage(sky, sky_level) + sky += roman.thermal_backgrounds[galsim_filter_name] * pixscalefrac ** 2 + else: + sky = sky * pixscalefrac ** 2 / etomjysr + # convert to electrons / s / output pixel + + # Flux in AB mags to electrons + maggytoes = romanisim.bandpass.get_abflux(filter_name) # Set effective read noise if effreadnoise is None: - effreadnoise = efftimes / romanisim.parameters.read_time + readnoise = np.median(romanisim.parameters.reference_data['readnoise']) + gain = np.median(romanisim.parameters.reference_data['gain']) + effreadnoise = ( + np.sqrt(2) * readnoise * gain) + # sqrt(2) from subtracting one read from another + effreadnoise /= (np.median(efftimes * pixscalefrac ** 2) / nexposures) + # divided by the typical exposure length + # the efftimes are multiplied by the square of the pixscalefrac + # to reflect the fact that if pixscalefrac < 1, each output pixel + # sees less of the total exposure time than the input pixels. + effreadnoise /= np.sqrt(nexposures) + # averaging down like the sqrt of the number of exposures + # note that we are ignoring all of the individual reads, which also + # contribute to reducing the effective read noise. Pass --effreadnoise + # if you want to do better than this! + effreadnoise = effreadnoise.to(u.electron).value * etomjysr + # converting to MJy/sr units + else: + effreadnoise = 0 + + chromatic = False + if (len(catalog) > 0 and not isinstance(catalog, astropy.table.Table) + and catalog[0].profile.spectral): + chromatic = True + + if psf is None: + if (pixscalefrac > 1) or (pixscalefrac < 0): + raise ValueError('weird pixscale!') + psf = l3_psf(filter_name, pixscalefrac, webbpsf=webbpsf, + chromatic=chromatic) # Simulate mosaic cps - mosaic, simcatobj = simulate_cps( - image, meta, efftimes, objlist=catalog, psf=psf, + mosaic, extras = simulate_cps( + image, filter_name, efftimes, objlist=catalog, psf=psf, sky=sky, - effreadnoise=effreadnoise, cps_conv=cps_conv, coords_unit=coords_unit, - bandpass=bandpass, - wcs=wcs, rng=rng, seed=seed, unit_factor=unit_factor, + effreadnoise=effreadnoise, bandpass=bandpass, + rng=rng, seed=seed, + maggytoes=maggytoes, etomjysr=etomjysr, **kwargs) - # Set poisson variance - if "var_poisson" in simcatobj: - var_poisson = simcatobj["var_poisson"] - else: - var_poisson = None - - # Set read noise variance - if "var_readnoise" in simcatobj: - var_readnoise = simcatobj["var_readnoise"].value - else: - var_readnoise = None - # Create Mosaic Model - mosaic_mdl = make_l3(mosaic, metadata, efftimes, unit_factor=unit_factor, - var_poisson=var_poisson, var_readnoise=var_readnoise) - - # Add simulation artifacts - extras = {} - extras['simcatobj'] = simcatobj - extras['wcs'] = mosaic.wcs + var_poisson = extras.pop('var_poisson') + var_rnoise = extras.pop('var_rnoise') + context = np.ones((1,) + mosaic.array.shape, dtype=np.uint32) + mosaic_mdl = make_l3(mosaic, meta, efftimes, var_poisson=var_poisson, + var_rnoise=var_rnoise, context=context) log.info('Simulation complete.') # return mosaic, extras return mosaic_mdl, extras -def simulate_cps(image, metadata, efftimes, objlist=None, psf=None, - wcs=None, xpos=None, ypos=None, coord=None, sky=0, - effreadnoise=None, cps_conv=1, unit_factor=1.0 * (u.MJy / u.sr), - bandpass=None, coords_unit='rad', - rng=None, seed=None, ignore_distant_sources=10,): +def get_pixscalefrac(wcs, shape): + """Compute pixel scale from WCS, scaled to nominal units of 0.11". + + Computes the difference in arcseconds between the central pixels. + Assumes square pixels. Scales this so that a nominal WCS with 0.11" + pixels gives a value of 1. + + Parameters + ---------- + wcs : galsim.WCS + WCS to get pixel scale for + shape : tuple + image shape + """ + cenpix = shape[0] // 2, shape[1] // 2 + pos1 = wcs.toWorld(galsim.PositionD(cenpix[0], cenpix[1])) + pos2 = wcs.toWorld(galsim.PositionD(cenpix[0], (cenpix[1] + 1))) + pixscale = pos1.distanceTo(pos2).deg * 60 * 60 # arcsec + pixscale /= romanisim.parameters.pixel_scale + if (pixscale > 1 and pixscale < 1.05): + pixscale = 1 + log.info('Setting pixscale to match default.') + return pixscale + + +def simulate_cps(image, filter_name, efftimes, objlist=None, psf=None, + xpos=None, ypos=None, coord=None, sky=0, bandpass=None, + effreadnoise=None, maggytoes=None, etomjysr=None, + rng=None, seed=None, ignore_distant_sources=10,): """Simulate average MegaJankies per steradian in a single SCA. Parameters ---------- image : galsim.Image Image onto which other effects should be added, with associated WCS. - metadata : dict - Metadata structure for Roman asdf file. + filter_name : str + filter to simulate efftimes : np.ndarray or float - Effective exposure time of reach pixel in mosaic. + Time Roman spent observing each part of the sky. If an array, shape must match shape parameter. objlist : list[CatalogObject], Table, or None Sources to render psf : galsim.Profile PSF to use when rendering sources - wcs : galsim.GSFitsWCS - WCS corresponding to image xpos, ypos : array_like (float) x, y positions of each source in objlist coord : array_like (float) - ra, dec positions of each source in objlist + ra, dec positions of each source in objlist (deg) sky : float or array_like - Image or constant with the counts / pix / sec from sky. - effreadnoise : float - Effective read noise for mosaic. - cps_conv : float - Factor to convert data to cps - unit_factor: float - Factor to convert counts data to MJy / sr + Image or constant with the electron / pix / sec from sky. bandpass : galsim.Bandpass - bandpass to use for rendering chromatic objects - coords_unit : string - units of coords + bandpass being used. Only used for chromatic objects + effreadnoise : float + Effective read noise for mosaic (MJy / sr) + maggytoes: float + Factor to convert electrons to MJy / sr; one maggy makes + this many e/s. + etomjysr : float + Factor to convert electron to MJy/sr; one e/s/pix corresponds + to this MJy/sr. rng : galsim.BaseDeviate random number generator seed : int @@ -344,15 +501,10 @@ def simulate_cps(image, metadata, efftimes, objlist=None, psf=None, ignore_distant_sources : int Ignore sources more than this distance off image. - Returns - ------- - objinfo : np.ndarray - Information on position and flux of each rendered source. - Returns ------- image : galsim.Image - Idealized image of scene as seen by Roman + Idealized image of scene as seen by Roman (MJy / sr) extras : dict catalog of simulated objects in image, noise, and misc. debug """ @@ -364,6 +516,12 @@ def simulate_cps(image, metadata, efftimes, objlist=None, psf=None, if rng is None: rng = galsim.UniformDeviate(seed) + if etomjysr is None: + etomjysr = romanisim.bandpass.etomjysr(filter_name) + + if maggytoes is None: + maggytoes = romanisim.bandpass.get_abflux(filter_name) + # Dictionary to hold simulation artifacts extras = {} @@ -376,180 +534,128 @@ def simulate_cps(image, metadata, efftimes, objlist=None, psf=None, objlist = [] if len(objlist) > 0 and xpos is None: if isinstance(objlist, table.Table): - coord = np.array([[o['ra'], np.radians(o['dec'])] for o in objlist]) + coord = np.array([[o['ra'], o['dec']] for o in objlist]) else: - if (coords_unit == 'rad'): - coord = np.array([[o.sky_pos.ra.rad, o.sky_pos.dec.rad] - for o in objlist]) - else: - coord = np.array([[o.sky_pos.ra.deg, o.sky_pos.dec.deg] - 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) - - chromatic = False - if len(objlist) > 0 and objlist[0].profile.spectral: - chromatic = True + coord = np.array([[o.sky_pos.ra.deg, o.sky_pos.dec.deg] + for o in objlist]) + xpos, ypos = image.wcs.radecToxy(coord[:, 0], coord[:, 1], 'deg') # Check for objects outside the image boundary (+ consideration) if len(objlist) > 0: - keep = romanisim.image.in_bounds(xpos, ypos, image.bounds, ignore_distant_sources) + xpos = np.array(xpos) + ypos = np.array(ypos) + keep = romanisim.image.in_bounds(xpos, ypos, image.bounds, + ignore_distant_sources) - objlist = np.array(objlist)[keep].tolist() + if isinstance(objlist, astropy.table.Table): + objlist = objlist[keep] + else: + objlist = [o for (o, k) in zip(objlist, keep) if k] xpos = xpos[keep] ypos = ypos[keep] + + if len(objlist) > 0: # Pixelized object locations - xpos_idx = [round(x) for x in xpos] - ypos_idx = [round(y) for y in ypos] + xpos_idx = np.round(xpos).astype('i4') + ypos_idx = np.round(ypos).astype('i4') offedge = romanisim.image.in_bounds(xpos, ypos, image.bounds, 0) # Set exposure time per source if isinstance(efftimes, np.ndarray): - src_exptimes = [efftimes[y, x] for x, y in zip(xpos_idx, ypos_idx)] + src_exptimes = [ + efftimes[y, x] if onframe else -1 + for x, y, onframe in zip(xpos_idx, ypos_idx, ~offedge)] + avg_exptime = np.average(efftimes[efftimes > 0]) else: src_exptimes = [efftimes] * len(xpos) + avg_exptime = efftimes src_exptimes = np.array(src_exptimes) - - # Set the average exposure time to objects lacking one - if True in offedge: - avg_exptime = np.average(src_exptimes) - src_exptimes[offedge] = avg_exptime - - else: - src_exptimes = [] - - # Generate GWCS compatible wcs - if wcs is None: - wcs = romanisim.wcs.get_mosaic_wcs(metadata, shape=image.array.shape, xpos=xpos, ypos=ypos, coord=coord) - image.wcs = wcs - - # Add objects to mosaic - if len(src_exptimes) > 0: - objinfokeep = add_objects_to_l3( - image, objlist, src_exptimes, wcs=wcs, xpos=xpos, ypos=ypos, - filter_name=metadata['basic']['optical_element'], bandpass=bandpass, psf=psf, rng=rng, - cps_conv=cps_conv, unit_factor=unit_factor) - - # Add object info artifacts - objinfo = {} - objinfo['array'] = np.zeros( - len(objlist), - dtype=[('x', 'f4'), ('y', 'f4'), ('counts', 'f4'), ('time', 'f4')]) - if len(objlist) > 0: + src_exptimes[src_exptimes == -1] = avg_exptime + + if isinstance(objlist, astropy.table.Table): + objlist = romanisim.catalog.table_to_catalog(objlist, [filter_name]) + + # Add objects to mosaic + chromatic = objlist[0].profile.spectral + maggytoes0 = maggytoes if not chromatic else 1 + objinfo0 = add_objects_to_l3( + image, objlist, src_exptimes, xpos=xpos, ypos=ypos, + filter_name=filter_name, psf=psf, bandpass=bandpass, rng=rng, + maggytoes=maggytoes0, etomjysr=etomjysr) + objinfo = np.zeros( + len(objlist), + dtype=[('x', 'f4'), ('y', 'f4'), ('counts', 'f4'), ('time', 'f4')]) objinfo['x'] = xpos objinfo['y'] = ypos - objinfo['counts'] = objinfokeep['counts'] - objinfo['time'] = objinfokeep['time'] - else: - objinfo['counts'] = np.array([]) - objinfo['time'] = np.array([]) - extras['objinfo'] = objinfo - - # Noise - im_no_noise = image.array.copy() - - # add Poisson noise if we made a noiseless, not-photon-shooting - # image. - poisson_noise = galsim.PoissonNoise(rng) - if not chromatic: - # This works in ADU, so need to convert back to counts first.. add this, then convert back - image *= efftimes / unit_factor.value - image.addNoise(poisson_noise) - image /= efftimes / unit_factor.value + objinfo['counts'] = objinfo0['counts'] + objinfo['time'] = objinfo0['time'] + extras['objinfo'] = objinfo if sky is not None: - if isinstance(sky, (galsim.Image)): - sky_arr = sky.array - elif not isinstance(sky, (np.ndarray)): - sky_arr = sky * np.ones(image.array.shape, dtype=float) - else: - sky_arr = sky - + # in e / s / output pixel + poisson_noise = galsim.PoissonNoise(rng) workim = image * 0 - workim += sky - workim *= efftimes + workim += sky * efftimes workim.addNoise(poisson_noise) - - image += (workim * unit_factor.value / efftimes) - - else: - sky_arr = np.zeros(image.array.shape, dtype=float) - - extras['var_poisson'] = np.abs(image.array - (im_no_noise + sky_arr)) + image += (workim * etomjysr / efftimes) # Add readnoise if effreadnoise is not None: - # This works in ADU, so need to convert back to counts first.. add this, then convert back + # in MJy / sr readnoise = np.zeros(image.array.shape, dtype='f4') rn_rng = galsim.GaussianDeviate(seed) rn_rng.generate(readnoise) readnoise = readnoise * effreadnoise - readnoise /= np.sqrt(efftimes) + image += readnoise + else: + effreadnoise = 0 + extras['var_rnoise'] = effreadnoise ** 2 - image += readnoise * unit_factor.value / efftimes - extras['readnoise'] = readnoise * unit_factor / efftimes + var_poisson_factor = (efftimes / etomjysr) * etomjysr ** 2 / efftimes ** 2 + # goofy game with etomjysr: image * (efftimes / etomjysr) + # -> number of photons entering each pixel + # then we interpret this as a variance (since mean = variance for a + # Poisson distribution), and convert the variance image + # to the final units with two factors of etomjysr and the effective time - extras['var_readnoise'] = (rn_rng.sigma * (effreadnoise / np.sqrt(efftimes)) * (unit_factor / efftimes))**2 + extras['var_poisson'] = ( + np.clip(image.array, 0, np.inf) * var_poisson_factor) # Return image and artifacts return image, extras -def make_l3(image, metadata, efftimes, rng=None, seed=None, var_poisson=None, - var_flat=None, var_readnoise=None, context=None, - unit_factor=(1.0 * u.MJy / u.sr)): +def make_l3(image, metadata, efftimes, var_poisson=None, + var_flat=None, var_rnoise=None, context=None): """ Create and populate MosaicModel of image and noises. Parameters ---------- image : galsim.Image - Image containing mosaic data. + Image containing mosaic data (MJy / sr) metadata : dict Metadata structure for Roman asdf file. efftimes : np.ndarray or float - Effective exposure time of reach pixel in mosaic. + Time Roman spent observing each part of the sky. If an array, shape must match shape parameter. - rng : galsim.BaseDeviate - Random number generator. - seed : int - Seed for random number generator. var_poisson : np.ndarray Poisson variance for each pixel var_flat : np.ndarray Flat variance for each pixel - var_readnoise : np.ndarray + var_rnoise : np.ndarray Read Noise variance for each pixel context : np.ndarray File number(s) for each pixel - unit_factor: float - Factor to convert counts data to MJy / sr - - Returns - ------- - objinfo : np.ndarray - Information on position and flux of each rendered source. Returns ------- - image : rdm.MosaicModel - Mosaic model object containing the data, metadata, variances, weight, and context. + image : roman_datamodels.datamodels.MosaicModel + Mosaic datamodel """ # Create mosaic data object - 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) + mosaic = image.array.copy() * u.MJy / u.sr # Ensure that effective times are an array if isinstance(efftimes, np.ndarray): @@ -558,43 +664,111 @@ def make_l3(image, metadata, efftimes, rng=None, seed=None, var_poisson=None, efftimes_arr = efftimes * np.ones(mosaic.shape, dtype=np.float32) # Set mosaic to be a mosaic node - mosaic_node = maker_utils.mk_level3_mosaic(shape=mosaic.shape, meta=metadata) + mosaic_node = maker_utils.mk_level3_mosaic( + shape=mosaic.shape, meta=metadata) + mosaic_node.meta.wcs = metadata['wcs'] # Set data - mosaic_node.data = mosaic * unit_factor.unit - - # Poisson noise variance - if var_poisson is None: - mosaic_node.var_poisson = (unit_factor**2 / efftimes_arr**2).astype(np.float32) - else: - mosaic_node.var_poisson = u.Quantity(var_poisson.astype(np.float32), unit_factor.unit ** 2) - - # Read noise variance - if var_readnoise is None: - mosaic_node.var_rnoise = u.Quantity(np.zeros(mosaic.shape, dtype=np.float32), unit_factor.unit ** 2) - else: - mosaic_node.var_rnoise = u.Quantity((var_readnoise * np.ones(mosaic.shape)).astype(np.float32), unit_factor.unit ** 2) - - # Flat variance - if var_flat is None: - mosaic_node.var_flat = u.Quantity(np.zeros(mosaic.shape, dtype=np.float32), unit_factor.unit ** 2) - else: - mosaic_node.var_flat = var_flat - - # Total error - mosaic_node.err = np.sqrt(mosaic_node.var_rnoise + mosaic_node.var_flat + mosaic_node.var_poisson).astype(np.float32) + mosaic_node.data = mosaic + + var_poisson = 0 if var_poisson is None else var_poisson + var_rnoise = 0 if var_rnoise is None else var_rnoise + var_flat = 0 if var_flat is None else var_flat + context = (np.ones((1,) + mosaic.shape, dtype=np.uint32) + if context is None else context) + + mosaic_node.var_poisson[...] = var_poisson * mosaic.unit ** 2 + mosaic_node.var_rnoise[...] = var_rnoise * mosaic.unit ** 2 + mosaic_node.var_flat[...] = var_flat * mosaic.unit ** 2 + mosaic_node.err[...] = np.sqrt( + var_poisson + var_rnoise + var_flat) * mosaic.unit + mosaic_node.context = context + # Weight # Use exptime weight + # could scale this down by pixfrac ** 2 mosaic_node.weight = efftimes_arr.astype(np.float32) - # Context - # Binary index of images that contribute to the pixel - # Defined by geometry.. if not, all non zero = 1. - if context is None: - mosaic_node.context = np.ones((1,) + mosaic.shape, dtype=np.uint32) - else: - mosaic_node.context = context - - # Return mosaic return mosaic_node + + +def add_more_metadata(metadata, efftimes, filter_name, wcs, shape, nexposures): + """Fill in the L3 metadata for simulations. + + Updates the 'metadata' array in place. Touches a number of fields in + metadata.basic, metadata.photometry, metadata.resample, metadata.wcsinfo + and the metadata root. + + Parameters + ---------- + metadata + metadata to update + efftimes : np.ndarray + exposure time on each pixel + filter_name : str + name of filter + wcs : gwcs.wcs.WCS + WCS for mosaic + shape : tuple[2] + shape of mosaic + nexposures : int + number of exposures contributing to mosaic + """ + maxtime = np.max(efftimes) + meantime = metadata['basic']['time_mean_mjd'] + meanexptime = (efftimes if np.isscalar(efftimes) else + np.mean(efftimes[efftimes > 0])) + # guesses at the first and last times; do not really make sense + # for this kind of simulation + metadata['basic']['time_first_mjd'] = meantime - maxtime / 24 / 60 / 60 / 2 + metadata['basic']['time_last_mjd'] = meantime + maxtime / 24 / 60 / 60 / 2 + metadata['basic']['max_exposure_time'] = maxtime + metadata['basic']['mean_exposure_time'] = meanexptime + for step in ['flux', 'outlier_detection', 'skymatch', 'resample']: + metadata['cal_step'][step] = 'COMPLETE' + metadata['basic']['individual_image_meta'] = None + metadata['model_type'] = 'WfiMosaic' + metadata['photometry']['conversion_microjanskys'] = ( + (1e12 * (u.rad / u.arcsec) ** 2).to(u.dimensionless_unscaled) * + u.uJy / u.arcsec ** 2) + metadata['photometry']['conversion_megajanskys'] = 1 * u.MJy / u.sr + + cenx, ceny = ((shape[1] - 1) / 2, (shape[0] - 1) / 2) + c1 = wcs.pixel_to_world(cenx, ceny) + c2 = wcs.pixel_to_world(cenx + 1, ceny) + pscale = c1.separation(c2) + + metadata['photometry']['pixelarea_steradians'] = (pscale ** 2).to(u.sr) + metadata['photometry']['pixelarea_arcsecsq'] = ( + pscale.to(u.arcsec) ** 2) + metadata['photometry']['conversion_microjanskys_uncertainty'] = ( + 0 * u.uJy / u.arcsec ** 2) + metadata['photometry']['conversion_megajanskys_uncertainty'] = ( + 0 * u.MJy / u.sr) + metadata['resample']['pixel_scale_ratio'] = ( + pscale.to(u.arcsec).value / romanisim.parameters.pixel_scale) + metadata['resample']['pixfrac'] = 0 + # our simulations sort of imply idealized 0 droplet size + metadata['resample']['pointings'] = nexposures + metadata['resample']['product_exposure_time'] = ( + metadata['basic']['max_exposure_time']) + xref, yref = wcs.world_to_pixel( + metadata['wcsinfo']['ra_ref'], metadata['wcsinfo']['dec_ref']) + metadata['wcsinfo']['x_ref'] = xref + metadata['wcsinfo']['y_ref'] = yref + metadata['wcsinfo']['rotation_matrix'] = [[1, 0], [0, 1]] + metadata['wcsinfo']['pixel_scale'] = pscale.to(u.arcsec).value + metadata['wcsinfo']['pixel_scale_local'] = metadata['wcsinfo']['pixel_scale'] + metadata['wcsinfo']['s_region'] = romanisim.wcs.create_s_region(wcs, shape) + metadata['wcsinfo']['pixel_shape'] = shape + metadata['wcsinfo']['ra_center'] = c1.ra.to(u.degree).value + metadata['wcsinfo']['dec_center'] = c1.dec.to(u.degree).value + xcorn, ycorn = [[0, shape[1] - 1, shape[1] - 1, 0], + [0, 0, shape[0] - 1, shape[0] - 1]] + ccorn = wcs.pixel_to_world(xcorn, ycorn) + for i, corn in enumerate(ccorn): + metadata['wcsinfo']['ra_corn{i+1}'] = corn.ra.to(u.degree).value + metadata['wcsinfo']['dec_corn{i+1}'] = corn.dec.to(u.degree).value + metadata['wcsinfo']['orientat_local'] = 0 + metadata['wcsinfo']['orientat'] = 0 diff --git a/romanisim/parameters.py b/romanisim/parameters.py index 6619ed3b..5728a86c 100644 --- a/romanisim/parameters.py +++ b/romanisim/parameters.py @@ -5,8 +5,6 @@ from astropy.time import Time from astropy import units as u -pixel_scale = 0.11 - read_pattern = {1: [[1 + x for x in range(8)], [9 + x for x in range(8)], [17 + x for x in range(8)], @@ -67,13 +65,6 @@ }, 'wcsinfo': {'ra_ref': 270.0, 'dec_ref': 66.0, - 'v2_ref': 0, - 'v3_ref': 0, - 'roll_ref': 0, - 'vparity': -1, - 'v3yangle': -60.0, - # I don't know what vparity and v3yangle should really be, - # but they are always -1 and -60 in existing files. }, } @@ -86,17 +77,6 @@ "linearity": None, "readnoise": 5.0 * u.DN, "saturation": 55000 * u.DN, - # Taken from roman_wfi_photom_0046.asdf - "photom": {"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, - } - } } nborder = 4 # number of border pixels used for reference pixels. @@ -172,3 +152,6 @@ # angle of V3 relative to +Y V3IdlYAngle = -60 + +# fiducial WFI pixel scale in arcseconds +pixel_scale = 0.11 diff --git a/romanisim/psf.py b/romanisim/psf.py index a669751f..db104a85 100644 --- a/romanisim/psf.py +++ b/romanisim/psf.py @@ -24,7 +24,7 @@ def make_one_psf(sca, filter_name, wcs=None, webbpsf=True, pix=None, - chromatic=False, oversample=4, **kw): + chromatic=False, oversample=4, extra_convolution=None, **kw): """Make a PSF profile for Roman at a specific detector location. Can construct both PSFs using galsim's built-in galsim.roman.roman_psfs @@ -43,6 +43,8 @@ def make_one_psf(sca, filter_name, wcs=None, webbpsf=True, pix=None, pixel location of PSF on focal plane oversample : int oversampling with which to sample WebbPSF PSF + extra_convolution : galsim.gsobject.GSObject or None + Additional convolution to add to PSF **kw : dict Additional keywords passed to galsim.roman.getPSF or webbpsf.calc_psf, depending on whether webbpsf is set. @@ -67,8 +69,11 @@ def make_one_psf(sca, filter_name, wcs=None, webbpsf=True, pix=None, filter_name = None defaultkw.update(**kw) scapos = galsim.PositionD(*pix) if pix is not None else None - return roman.getPSF(sca, filter_name, wcs=wcs, SCA_pos=scapos, - wavelength=bandpass, **defaultkw) + res = roman.getPSF(sca, filter_name, wcs=wcs, SCA_pos=scapos, + wavelength=bandpass, **defaultkw) + if extra_convolution is not None: + res = galsim.Convolve(res, extra_convolution) + return res if chromatic: log.warning('romanisim does not yet support chromatic PSFs ' 'with webbpsf') @@ -115,11 +120,15 @@ def make_one_psf(sca, filter_name, wcs=None, webbpsf=True, pix=None, centroid = None intimg = galsim.InterpolatedImage( gimg, normalization='flux', use_true_center=True, offset=centroid) + + if extra_convolution is not None: + intimg = galsim.Convolve(intimg, extra_convolution) + return intimg def make_psf(sca, filter_name, wcs=None, webbpsf=True, pix=None, - chromatic=False, variable=False, **kw): + chromatic=False, variable=False, extra_convolution=None, **kw): """Make a PSF profile for Roman. Optionally supports spatially variable PSFs via interpolation between @@ -138,6 +147,8 @@ def make_psf(sca, filter_name, wcs=None, webbpsf=True, pix=None, pixel location of PSF on focal plane variable : bool True if a variable PSF object is desired + extra_convolution : galsim.gsobject.GSObject or None + Additional convolution to add to PSF profiles **kw : dict Additional keywords passed to make_one_psf @@ -149,7 +160,8 @@ def make_psf(sca, filter_name, wcs=None, webbpsf=True, pix=None, """ if not variable: return make_one_psf(sca, filter_name, wcs=wcs, webbpsf=webbpsf, - pix=pix, chromatic=chromatic, **kw) + pix=pix, chromatic=chromatic, + extra_convolution=extra_convolution, **kw) elif pix is not None: raise ValueError('cannot set both pix and variable') buf = 49 @@ -162,7 +174,8 @@ def make_psf(sca, filter_name, wcs=None, webbpsf=True, pix=None, psfs = dict() for corner, pix in corners.items(): psfs[corner] = make_one_psf(sca, filter_name, wcs=wcs, webbpsf=webbpsf, - pix=pix, chromatic=chromatic, **kw) + pix=pix, chromatic=chromatic, + extra_convolution=extra_convolution, **kw) return VariablePSF(corners, psfs) diff --git a/romanisim/tests/test_image.py b/romanisim/tests/test_image.py index 717e23c6..accc6929 100644 --- a/romanisim/tests/test_image.py +++ b/romanisim/tests/test_image.py @@ -30,7 +30,6 @@ from metrics_logger.decorators import metrics_logger from romanisim import log from roman_datamodels.stnode import WfiScienceRaw, WfiImage -from roman_datamodels import datamodels import romanisim.bandpass @@ -654,7 +653,6 @@ def test_inject_source_into_image(): meta, cat, usecrds=False, webbpsf=True, level=2, rng=rng, psf_keywords=dict(nlambda=1), crparam=None) - im = datamodels.ImageModel(im) # Create catalog with one source for injection xpos, ypos = 10, 10 diff --git a/romanisim/tests/test_l3.py b/romanisim/tests/test_l3.py index a1238330..0d828c0b 100644 --- a/romanisim/tests/test_l3.py +++ b/romanisim/tests/test_l3.py @@ -10,6 +10,7 @@ from romanisim import parameters, catalog, wcs, l3, psf, util, log from astropy import units as u from astropy import table +from astropy.stats import mad_std import asdf import pytest from metrics_logger.decorators import metrics_logger @@ -33,7 +34,8 @@ def test_inject_sources_into_mosaic(): metadata['basic']['optical_element'] = filter_name # Create WCS - twcs = wcs.get_mosaic_wcs(metadata, shape=(galsim.roman.n_pix, galsim.roman.n_pix)) + twcs = wcs.GWCS(wcs.get_mosaic_wcs( + metadata, shape=(galsim.roman.n_pix, galsim.roman.n_pix))) # Create initial Level 3 mosaic @@ -42,14 +44,18 @@ def test_inject_sources_into_mosaic(): # (total files contributed to each quadrant) # Create gaussian noise generators - g1 = galsim.GaussianDeviate(rng_seed, mean=1.0e-7, sigma=0.01e-7) - g2 = galsim.GaussianDeviate(rng_seed, mean=1.0e-7, sigma=0.02e-7) - g3 = galsim.GaussianDeviate(rng_seed, mean=1.0e-7, sigma=0.05e-7) - g4 = galsim.GaussianDeviate(rng_seed, mean=1.0e-7, sigma=0.1e-7) + # sky should generate ~0.2 electron / s / pix. + # MJy / sr has similar magnitude to electron / s (i.e., within a factor + # of several), so just use 0.2 here. + meanflux = 0.2 + g1 = galsim.GaussianDeviate(rng_seed, mean=meanflux, sigma=0.01 * meanflux) + g2 = galsim.GaussianDeviate(rng_seed, mean=meanflux, sigma=0.02 * meanflux) + g3 = galsim.GaussianDeviate(rng_seed, mean=meanflux, sigma=0.05 * meanflux) + g4 = galsim.GaussianDeviate(rng_seed, mean=meanflux, sigma=0.10 * meanflux) # 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 + l3_mos['meta']['wcs'] = twcs._wcs # Update metadata in the l3 model for key in metadata.keys(): @@ -57,12 +63,10 @@ def test_inject_sources_into_mosaic(): l3_mos.meta[key].update(metadata[key]) # Obtain unit conversion factors - # Need to convert from counts / pixel to MJy / sr - # Flux to counts + # maggies to counts (large number) cps_conv = romanisim.bandpass.get_abflux(filter_name) - # Unit factor - unit_factor = ((3631 * u.Jy) / (romanisim.bandpass.get_abflux(filter_name) * 10e6 - * parameters.reference_data['photom']["pixelareasr"][filter_name])).to(u.MJy / u.sr) + # electrons to mjysr (roughly order unity in scale) + unit_factor = romanisim.bandpass.etomjysr(filter_name) # Populate the mosaic data array with gaussian noise from generators g1.generate(l3_mos.data.value[0:100, 0:100]) @@ -71,13 +75,13 @@ def test_inject_sources_into_mosaic(): 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 + l3_mos.var_poisson.value[0:100, 0:100] = 0.01**2 * meanflux**2 + l3_mos.var_poisson.value[0:100, 100:200] = 0.02**2 * meanflux**2 + l3_mos.var_poisson.value[100:200, 0:100] = 0.05**2 * meanflux**2 + l3_mos.var_poisson.value[100:200, 100:200] = 0.1**2 * meanflux**2 # Create normalized psf source catalog (same source in each quadrant) - mag_flux = 1e-10 + mag_flux = 1e-9 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], filter_name: 4 * [mag_flux]} sc_table = table.Table(sc_dict) @@ -86,19 +90,9 @@ def test_inject_sources_into_mosaic(): 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[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) - - 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, [filter_name]) + ra, dec = twcs._radec(np.array(xpos_idx), np.array(ypos_idx)) + sc_table['ra'] = np.degrees(ra) + sc_table['dec'] = np.degrees(dec) # Copy original Mosaic before adding sources as sources are added in place l3_mos_orig = l3_mos.copy() @@ -106,15 +100,8 @@ 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, cps_conv=cps_conv, 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, - 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.inject_sources_into_l3( + l3_mos, sc_table, seed=rng_seed) # Ensure that every data pixel value has increased or # remained the same with the new sources injected @@ -128,8 +115,8 @@ def test_inject_sources_into_mosaic(): assert np.all(l3_mos.var_poisson.value[~close_mask] > l3_mos_orig.var_poisson.value[~close_mask]) # Ensure total added flux matches expected added flux - total_rec_flux = np.sum(l3_mos.data - l3_mos_orig.data) / unit_factor - total_theo_flux = 4 * mag_flux * cps_conv + total_rec_flux = np.sum(l3_mos.data - l3_mos_orig.data) # MJy / sr + total_theo_flux = 4 * mag_flux * cps_conv * unit_factor * u.MJy / u.sr assert np.isclose(total_rec_flux, total_theo_flux, rtol=4e-02) # Create log entry and artifacts @@ -169,13 +156,10 @@ def test_sim_mosaic(): cat = catalog.make_dummy_table_catalog(cen, radius=0.02, nobj=100, seed=rng_seed) # Make the first 10 bright for tests cat[filter_name][0:10] *= 1e4 - source_cat = catalog.table_to_catalog(cat, [filter_name]) # Create bounds from the object list twcs = romanisim.wcs.get_mosaic_wcs(metadata) - coords = np.array([[o.sky_pos.ra.rad, o.sky_pos.dec.rad] - for o in source_cat]) - allx, ally = twcs.radecToxy(coords[:, 0], coords[:, 1], 'rad') + allx, ally = twcs.world_to_pixel(cat['ra'], cat['dec']) # Obtain the sample extremums xmin = min(allx) @@ -184,7 +168,7 @@ def test_sim_mosaic(): ymax = max(ally) # Obtain WCS center - xcen, ycen = twcs.radecToxy(twcs.center.ra, twcs.center.dec, 'rad') + xcen, ycen = twcs.world_to_pixel(ra_ref, dec_ref) # Determine maximum extremums from WCS center xdiff = max([math.ceil(xmax - xcen), math.ceil(xcen - xmin)]) + 1 @@ -198,21 +182,22 @@ def test_sim_mosaic(): # Simulate mosaic mosaic, extras = l3.simulate(context.shape[1:], moswcs, exptimes[0], - filter_name, source_cat, metadata=metadata, seed=rng_seed) + filter_name, cat, metadata=metadata, + seed=rng_seed) + + # Did all sources get simulated? + assert len(extras['objinfo']) == len(cat) # Ensure center pixel of bright objects is bright - x_all, y_all = moswcs.radecToxy(coords[:10, 0], coords[:10, 1], 'rad') + x_all, y_all = moswcs.world_to_pixel(cat['ra'][:10], cat['dec'][:10]) for x, y in zip(x_all, y_all): x = int(x) y = int(y) assert mosaic.data.value[y, x] > (np.median(mosaic.data.value) * 5) # Did we get all the flux? - # Convert to CPS for comparison - # Unit factor - unit_factor = ((3631 * u.Jy) / (romanisim.bandpass.get_abflux(filter_name) * 10e6 - * parameters.reference_data['photom']["pixelareasr"][filter_name])).to(u.MJy / u.sr) - totflux = np.sum(mosaic.data.value - np.median(mosaic.data.value)) / unit_factor.value + etomjysr = romanisim.bandpass.etomjysr(filter_name) + totflux = np.sum(mosaic.data.value - np.median(mosaic.data.value)) / etomjysr # Flux to counts cps_conv = romanisim.bandpass.get_abflux(filter_name) @@ -221,6 +206,14 @@ def test_sim_mosaic(): # Ensure that the measured flux is close to the expected flux assert np.abs(np.log(expectedflux) / np.log(totflux) - 1) < 0.1 + # Is the noise about right? + assert np.abs( + mad_std(mosaic.data.value) / np.median(mosaic.err.value) - 1) < 0.5 + # note large ~50% error bar there; value in initial test run is 0.25 + # a substantial number of source pixels have flux, so the simple medians + # and mads aren't terribly right. + # if I repeat this after only including the first source I get 1.004. + # Add log entries and artifacts log.info('DMS219 successfully created mosaic file with sources rendered ' 'at correct locations with matching flux and added noise.') @@ -246,13 +239,16 @@ def set_up_image_rendering_things(): chromatic=True) bandpass = roman.getBandpasses(AB_zeropoint=True)['H158'] counts = 1000 + maggiestoe = romanisim.bandpass.get_abflux(filter_name) fluxdict = {filter_name: counts} + fluxdictgray = {filter_name: counts / maggiestoe} # Sample catalogs graycatalog = [ - catalog.CatalogObject(None, galsim.DeltaFunction(), deepcopy(fluxdict)), + catalog.CatalogObject(None, galsim.DeltaFunction(), + deepcopy(fluxdictgray)), catalog.CatalogObject(None, galsim.Sersic(1, half_light_radius=0.2), - deepcopy(fluxdict)) + deepcopy(fluxdictgray)) ] vega_sed = galsim.SED('vega.txt', 'nm', 'flambda') vega_sed = vega_sed.withFlux(counts, bandpass) @@ -264,7 +260,7 @@ def set_up_image_rendering_things(): tabcat = table.Table() tabcat['ra'] = [270.0] tabcat['dec'] = [66.0] - tabcat[filter_name] = counts + tabcat[filter_name] = counts / maggiestoe tabcat['type'] = 'PSF' tabcat['n'] = -1 tabcat['half_light_radius'] = -1 @@ -315,35 +311,43 @@ def test_simulate_vs_cps(): im = imdict['im'].copy() im.array[:] = 0 - # Set WCS - twcs = romanisim.wcs.get_mosaic_wcs(metadata, shape=(roman.n_pix, roman.n_pix)) + maggytoes = romanisim.bandpass.get_abflux(filter_name) + etomjysr = romanisim.bandpass.etomjysr(filter_name) + + twcs = wcs.get_mosaic_wcs(meta, shape=im.array.shape) + im.wcs = wcs.GWCS(twcs) # Create chromatic data in simulate_cps + im1 = im.copy() - im1, extras1 = l3.simulate_cps(im1, metadata, exptime, objlist=chromcat, - xpos=[50] * len(chromcat), ypos=[50] * len(chromcat), - bandpass=imdict['bandpass'], seed=rng_seed, - ignore_distant_sources=100) + im1, extras1 = l3.simulate_cps(im1, filter_name, exptime, objlist=chromcat, + bandpass=imdict['bandpass'], + psf=imdict['impsfchromatic'], + seed=rng_seed, + ignore_distant_sources=100, + maggytoes=maggytoes, etomjysr=etomjysr) # Create filter data in simulate_cps im2 = im.copy() - im2, extras2 = l3.simulate_cps(im2, metadata, exptime, objlist=graycat, - xpos=[50] * len(chromcat), ypos=[50] * len(chromcat), + im2, extras2 = l3.simulate_cps(im2, filter_name, exptime, objlist=graycat, psf=imdict['impsfgray'], seed=rng_seed, - ignore_distant_sources=100) + ignore_distant_sources=100, + maggytoes=maggytoes, etomjysr=etomjysr) # Create chromatic data in simulate im3, extras3 = l3.simulate((roman.n_pix, roman.n_pix), twcs, exptime, filter_name, chromcat, - bandpass=imdict['bandpass'], seed=rng_seed, - cps_conv=1, unit_factor=(1 * u.MJy / u.sr), + bandpass=imdict['bandpass'], + psf=imdict['impsfchromatic'], + seed=rng_seed, metadata=metadata, sky=0, - ignore_distant_sources=100, effreadnoise=0, + ignore_distant_sources=100, + effreadnoise=0, ) # Create filter data in simulate im4, extras4 = l3.simulate((roman.n_pix, roman.n_pix), twcs, exptime, filter_name, graycat, psf=imdict['impsfgray'], - cps_conv=1, unit_factor=(1 * u.MJy / u.sr), seed=rng_seed, + seed=rng_seed, metadata=metadata, sky=0, ignore_distant_sources=100, effreadnoise=0, @@ -378,7 +382,7 @@ def test_simulate_cps(): # Test empty image l3.simulate_cps( - im, metadata, exptime, objlist=[], psf=imdict['impsfgray'], + im, filter_name, exptime, objlist=[], psf=imdict['impsfgray'], sky=0) assert np.all(im.array == 0) # verify nothing in -> nothing out @@ -387,7 +391,8 @@ def test_simulate_cps(): skycountspersecond = 1 sky.array[:] = skycountspersecond im2 = im.copy() - l3.simulate_cps(im2, metadata, exptime, sky=sky, seed=rng_seed) + l3.simulate_cps(im2, filter_name, exptime, sky=sky, seed=rng_seed, + etomjysr=1) # verify adding the sky increases the counts assert np.all(im2.array >= im.array) # verify that the count rate is about right. @@ -415,28 +420,28 @@ def test_simulate_cps(): # render some objects im3 = im.copy() _, objinfo = l3.simulate_cps( - im3, metadata, exptime, objlist=imdict['graycatalog'], psf=imdict['impsfgray'], + im3, filter_name, exptime, objlist=imdict['graycatalog'], psf=imdict['impsfgray'], xpos=[50, 50], ypos=[50, 50], seed=rng_seed) assert np.sum(im3.array) > 0 # at least verify that we added some sources... - assert len(objinfo['objinfo']['array']) == 2 # two sources were added + assert len(objinfo['objinfo']) == 2 # two sources were added im4 = im.copy() _, objinfo = l3.simulate_cps( - im4, metadata, exptime, objlist=imdict['chromcatalog'], + im4, filter_name, exptime, objlist=imdict['chromcatalog'], xpos=[50, 50], ypos=[50, 50], seed=rng_seed, psf=imdict['impsfchromatic'], bandpass=imdict['bandpass']) assert np.sum(im4.array) > 0 # at least verify that we added some sources... - assert len(objinfo['objinfo']['array']) == 2 # two sources were added + assert len(objinfo['objinfo']) == 2 # two sources were added im5 = im.copy() _, objinfo = l3.simulate_cps( - im5, metadata, exptime, objlist=imdict['chromcatalog'], + im5, filter_name, exptime, objlist=imdict['chromcatalog'], psf=imdict['impsfchromatic'], xpos=[1000, 1000], seed=rng_seed, ypos=[1000, 1000]) - assert np.sum(objinfo['objinfo']['counts'] > 0) == 0 + assert 'objinfo' not in objinfo # these sources should be out of bounds @@ -474,8 +479,10 @@ def test_exptime_array(): # Set variable exposure time array exptime = np.ones((roman.n_pix, roman.n_pix)) - exptime[0:50, :] = 300 - exptime[50:, :] = 600 + basetime = 300 + expfactor = 2 + exptime[0:50, :] = basetime + exptime[50:, :] = basetime * expfactor # Set WCS twcs = romanisim.wcs.get_mosaic_wcs(metadata, shape=(roman.n_pix, roman.n_pix)) @@ -483,24 +490,104 @@ def test_exptime_array(): # Create chromatic data simulation im1, extras1 = l3.simulate((roman.n_pix, roman.n_pix), twcs, exptime, filter_name, chromcat, bandpass=imdict['bandpass'], seed=rng_seed, - cps_conv=1, unit_factor=(1 * u.MJy / u.sr), metadata=metadata, ignore_distant_sources=100, + effreadnoise=0, ) # Create filter data simulation im2, extras2 = l3.simulate((roman.n_pix, roman.n_pix), twcs, exptime, filter_name, graycat, psf=imdict['impsfgray'], - cps_conv=1, unit_factor=(1 * u.MJy / u.sr), seed=rng_seed, metadata=metadata, ignore_distant_sources=100, + effreadnoise=0, ) # Ensure that the poisson variance scales with exposure time difference - assert np.isclose((np.mean(im1['var_poisson'][0:50, :].value) / np.mean(im1['var_poisson'][50:, :].value)), 2**0.5, rtol=0.1) - assert np.isclose((np.mean(im2['var_poisson'][0:50, :].value) / np.mean(im2['var_poisson'][50:, :].value)), 2**0.5, rtol=0.1) + assert np.isclose((np.median(im1['var_poisson'][0:50, :].value) / np.median(im1['var_poisson'][50:, :].value)), expfactor, rtol=0.02) + assert np.isclose((np.median(im2['var_poisson'][0:50, :].value) / np.median(im2['var_poisson'][50:, :].value)), expfactor, rtol=0.02) # Ensure that the data remains consistent across the exposure times - assert np.isclose(np.mean(im1['data'][0:50, :].value), np.mean(im1['data'][50:, :].value), rtol=0.2) - assert np.isclose(np.mean(im2['data'][0:50, :].value), np.mean(im2['data'][50:, :].value), rtol=0.2) + assert np.isclose(np.median(im1['data'][0:50, :].value), np.median(im1['data'][50:, :].value), rtol=0.02) + assert np.isclose(np.median(im2['data'][0:50, :].value), np.median(im2['data'][50:, :].value), rtol=0.02) -# TBD: Test of geometry construction + +def test_scaling(): + npix = 200 + imdict = set_up_image_rendering_things() + rng_seed = 1 + exptime = 400 + pscale = 0.1 + coord = SkyCoord(270 * u.deg, 66 * u.deg) + + # Set WCS + twcs1 = romanisim.wcs.create_tangent_plane_gwcs( + (npix / 2, npix / 2), pscale, coord) + twcs2 = romanisim.wcs.create_tangent_plane_gwcs( + (npix / 2, npix / 2), pscale / 2, coord) + + im1, extras1 = l3.simulate( + (npix, npix), twcs1, exptime, imdict['filter_name'], + imdict['tabcatalog'], seed=rng_seed, effreadnoise=0, + ) + + # half pixel scale + im2, extras2 = l3.simulate( + (npix * 2, npix * 2), twcs2, exptime, imdict['filter_name'], + imdict['tabcatalog'], seed=rng_seed, effreadnoise=0) + + # check that sky level doesn't depend on pixel scale (in calibrated units!) + skyfracdiff = np.median(im1.data.value) / np.median(im2.data.value) - 1 + log.info(f'skyfracdiff: {skyfracdiff:.3f}') + assert np.abs(skyfracdiff) < 0.1 + + # check that uncertainties match observed standard deviations + err1fracdiff = mad_std(im1.data.value) / np.median(im1.err.value) - 1 + err2fracdiff = mad_std(im2.data.value) / np.median(im2.err.value) - 1 + log.info(f'err1fracdiff: {err1fracdiff:.3f}, ' + f'err2fracdiff: {err2fracdiff:.3f}') + assert np.abs(err1fracdiff) < 0.1 + assert np.abs(err2fracdiff) < 0.1 + + # doubled exposure time + im3, extras3 = l3.simulate( + (npix, npix), twcs1, exptime * 10, imdict['filter_name'], + imdict['tabcatalog'], seed=rng_seed, effreadnoise=0) + + # check that sky level doesn't depend on exposure time (in calibrated units!) + sky3fracdiff = np.median(im1.data.value) / np.median(im3.data.value) - 1 + log.info(f'sky3fracdiff: {sky3fracdiff:.4f}') + assert np.abs(sky3fracdiff) < 0.1 + + # check that variances still work out + err3fracdiff = mad_std(im3.data.value) / np.median(im3.err.value) - 1 + log.info(f'err3fracdiff: {err3fracdiff:.3f}') + assert np.abs(err3fracdiff) < 0.1 + + # check that new variances are smaller than old ones by an appropriate factor + errfracdiff = ( + np.median(im1.err.value) / np.median(im3.err.value) - np.sqrt(10)) + log.info(f'err3 ratio diff from 1/sqrt(10) err1: {errfracdiff:0.3f}') + assert np.abs(errfracdiff) < 0.1 + + # check that fluxes match + # pixel scales are different by a factor of two. + fluxes = [] + for im, fac in zip((im1, im2, im3), (1, 2, 1)): + pix = im.meta.wcs.world_to_pixel( + imdict['tabcatalog']['ra'][0], imdict['tabcatalog']['dec'][0]) + pind = [int(x) for x in pix] + margin = 30 * fac + flux = np.sum(im.data.value[pind[1] - margin: pind[1] + margin, + pind[0] - margin: pind[0] + margin]) + fluxes.append(flux / fac ** 2) + # division by fac ** 2 accounts for different pixel scale + # i.e., we should be doing an integral here, and the pixels + # are a factor of 4 smaller in the second integral + # fluxes must match + for flux in fluxes[1:]: + fluxfracdiff = flux / fluxes[0] - 1 + log.info(f'fluxfracdiff: {fluxfracdiff:.5f}') + assert np.abs(fluxfracdiff) < 0.1 + + # these all match to a few percent; worst case in initial test run + # was err3fracdiff of 0.039. diff --git a/romanisim/wcs.py b/romanisim/wcs.py index 8c0f6abf..0ee540e5 100644 --- a/romanisim/wcs.py +++ b/romanisim/wcs.py @@ -459,23 +459,48 @@ def get_mosaic_wcs(mosaic, shape=None, xpos=None, ypos=None, coord=None): mosaic_node.data.shape[1]) shape = mosaic_node.data.shape - if (elem is None for elem in [xpos,ypos,coord]): - # Create a tangent plane WCS for the mosaic - # The affine parameters below should be reviewed and updated - affine = galsim.fitswcs.AffineTransform( - romanisim.parameters.pixel_scale, 0, 0, romanisim.parameters.pixel_scale, origin=galsim.PositionI(x=math.ceil(shape[1] / 2.0), y=math.ceil(shape[0] / 2.0)), - world_origin=galsim.PositionD(0, 0)) - wcs = galsim.fitswcs.TanWCS(affine, - util.celestialcoord(world_pos)) - else: - # Create GWCS compatible tangent plane WCS - header = {} - wcs = galsim.FittedSIPWCS(xpos, ypos, coord[:, 0], coord[:, 1], wcs_type='TAN', center=util.celestialcoord(world_pos)) - wcs._writeHeader(header, galsim.BoundsI(0, shape[0], 0, shape[1])) - # metadata['wcs'] = romanisim.wcs.wcs_from_fits_header(header) + xcen = math.ceil(shape[1] / 2.0) + ycen = math.ceil(shape[0] / 2.0) + wcs = create_tangent_plane_gwcs( + (xcen, ycen), romanisim.parameters.pixel_scale, world_pos) return wcs +def create_tangent_plane_gwcs(center, scale, center_coord): + """Create a tangent plane GWCS object. + + Parameters + ---------- + center : (float, float) + pixel coordinates of center + scale : float + pixel scale (arcsec) + center_coord : SkyCoord or CelestialCoord + coordinates of center pixel + + Returns + ------- + GWCS object corresponding to tangent plane projection + """ + + center_coord = util.skycoord(center_coord) + pixelshift = models.Shift(-center[0]) & models.Shift(-center[1]) + pixelscale = models.Scale(scale / 3600.) & models.Scale(scale / 3600.) + tangent_projection = models.Pix2Sky_TAN() + celestial_rotation = models.RotateNative2Celestial( + center_coord.ra.to(u.deg).value, + center_coord.dec.to(u.deg).value, 180.) + det2sky = pixelshift | pixelscale | tangent_projection | celestial_rotation + detector_frame = cf.Frame2D(name="detector", axes_names=("x", "y"), + unit=(u.pix, u.pix)) + sky_frame = cf.CelestialFrame( + reference_frame=astropy.coordinates.ICRS(), name='icrs', + unit=(u.deg, u.deg)) + wcsobj = gwcs.wcs.WCS([(detector_frame, det2sky), + (sky_frame, None)]) + return wcsobj + + def create_s_region(wcs, shape=None): """Create s_region string from wcs. diff --git a/scripts/romanisim-make-l3 b/scripts/romanisim-make-l3 new file mode 100755 index 00000000..703b35a3 --- /dev/null +++ b/scripts/romanisim-make-l3 @@ -0,0 +1,90 @@ +#!/usr/bin/env python + +import argparse +import yaml + +from astropy.coordinates import SkyCoord +from astropy import units as u +from astropy.table import Table +from astropy.time import Time +import galsim +import romanisim +from romanisim import log, wcs, persistence, parameters, l3, bandpass, util +from romanisim import ris_make_utils as ris +from copy import deepcopy +import math +import asdf +import os.path + + +if __name__ == '__main__': + parser = argparse.ArgumentParser( + description='Make an L3 image.', + epilog='EXAMPLE: %(prog)s output_image.asdf catalog.ecsv', + formatter_class=argparse.ArgumentDefaultsHelpFormatter) + parser.add_argument('filename', type=str, help='output image (asdf)') + parser.add_argument('catalog', type=str, help='input catalog (ecsv)') + parser.add_argument('--bandpass', type=str, help='bandpass to simulate', + default='F087') + parser.add_argument('--config', type=str, help='input parameter override file (yaml)', + default=None) + parser.add_argument('--date', type=str, default=None, + help=('UTC Date and Time of observation to simulate in ISOT format.')) + parser.add_argument('--radec', type=float, nargs=2, + help='ra and dec (deg)', default=None) + parser.add_argument('--npix', type=int, default=4000, + help='number of pixels across image') + parser.add_argument('--pixscalefrac', type=float, default=1.0, + help='pixel scale as fraction of original') + parser.add_argument('--exptime', type=float, default=100.0, + help='total exposure time on field; ' + 'roughly time per exposure times number of exposures') + parser.add_argument('--effreadnoise', type=float, default=None, + help='effective readnoise per pixel in MJy/sr. If ' + 'None, a pessimistic estimate is computed.') + parser.add_argument('--nexposures', type=float, default=1, + help='number of exposures on field. Used only ' + 'to compute the effective read noise.') + parser.add_argument('--rng_seed', type=int, default=None) + + args = parser.parse_args() + + log.info('Starting simulation...') + log.warning("romanisim is under active development. Its output has " + "not been formally validated; only limited testing has been " + "performed. For this reason, use of romanisim for " + "preparation of ROSES proposals is not advised. Other " + "packages like galsim's roman package or STIPS may better " + "serve such purposes.") + + pixscale = args.pixscalefrac * parameters.pixel_scale + midpoint = (args.npix - 1) / 2 + r, d = args.radec + center = util.celestialcoord(SkyCoord(ra=r * u.deg, dec=d * u.deg)) + twcs = wcs.create_tangent_plane_gwcs( + (midpoint, midpoint), pixscale, center) + cat = Table.read(args.catalog) + + metadata = dict() + if args.date is not None: + metadata['basic'] = dict(time_mean_mjd=Time(args.date, format='isot')) + metadata['filename'] = os.path.basename(args.filename) + + im, extras = l3.simulate( + (args.npix, args.npix), twcs, args.exptime, args.bandpass, + cat, effreadnoise=args.effreadnoise, nexposures=args.nexposures, + metadata=metadata) + + # Create metadata for simulation parameter + romanisimdict = deepcopy(vars(args)) + if 'filename' in romanisimdict: + romanisimdict['filename'] = str(romanisimdict['filename']) + romanisimdict.update(**extras) + romanisimdict['version'] = romanisim.__version__ + + af = asdf.AsdfFile() + af.tree = {'roman': im, 'romanisim': romanisimdict} + af.write_to(open(args.filename, 'wb')) + + +