7. Simulating Multiple Exposures](#mult_sims) \n",
"* [8. Extra time: Simulate deep field via galaxy catalog](#deep_field) \n",
"* [9. Calibrate the data](#calibrate_data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"\n",
"# 1. Getting Started\n",
"\n",
"\n",
"**Important:** \n",
"Before proceeding, ensure you have set the MIRAGE_DATA environment variable to point to the directory that contains the reference files associated with MIRAGE. \n",
"
\n",
"If you want JWST pipeline calibration reference files to be downloaded in a specific directory, you should also set the CRDS_DATA environment variable to point to that directory. This directory will also be used by the JWST calibration pipeline during data reduction.\n",
"
\n",
"You may also want to set the CRDS_SERVER_URL environment variable set to https://jwst-crds.stsci.edu. This is not strictly necessary, and Mirage will do it for you if you do not set it, but if you import the crds package, or any package that imports the crds package, you should set this environment variable first, in order to avoid an error.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import os"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# For use during JWebbinar\n",
"os.environ[\"MIRAGE_DATA\"] = \"/home/shared/mirage-data\"\n",
"os.environ[\"CRDS_DATA\"] = \"/home/jovyan/crds_cache\"\n",
"os.environ[\"CRDS_SERVER_URL\"] = \"https://jwst-crds.stsci.edu\"\n",
"\n",
"# Example when running notebook locally\n",
"#os.environ[\"MIRAGE_DATA\"] = \"/path/to/your/mirage_data\"\n",
"#os.environ[\"CRDS_DATA\"] = \"$HOME/crds_cache\"\n",
"#os.environ[\"CRDS_SERVER_URL\"] = \"https://jwst-crds.stsci.edu\""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# For examining outputs\n",
"from glob import glob\n",
"import numpy as np\n",
"from astropy.io import ascii, fits\n",
"from astropy.table import Table\n",
"import matplotlib.pyplot as plt\n",
"import urllib\n",
"import yaml\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# mirage imports\n",
"from mirage import imaging_simulator\n",
"from mirage.catalogs import create_catalog\n",
"from mirage.catalogs import catalog_generator\n",
"from mirage.seed_image import catalog_seed_image\n",
"from mirage.dark import dark_prep\n",
"from mirage.ramp_generator import obs_generator\n",
"from mirage.yaml import yaml_generator"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define the APT files that will be used for these simuations. As we will see throughout this notebook, the easiest way to create a consistent set of simulated exposures is to start with an APT program. Export the xml and pointing files from APT."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"xml_filename = 'sample_imaging.xml'\n",
"pointing_filename = 'sample_imaging.pointing'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Download the example xml and pointing files from APT\n",
"box_xml_file = 'https://stsci.box.com/shared/static/e8idi6u8yauvz1y8prwpe39d8e2lmigc.xml'\n",
"box_pointing_file = 'https://stsci.box.com/shared/static/izlyvihtzrefqzo5rn7zs2gmwg1lqb27.pointing'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"urllib.request.urlretrieve(box_xml_file, xml_filename)\n",
"urllib.request.urlretrieve(box_pointing_file, pointing_filename)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"# 2. Define convenience function for image display"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def show(array, title, min=0, max=1000):\n",
" plt.figure(figsize=(12, 12))\n",
" plt.imshow(array, clim=(min, max), origin='lower')\n",
" plt.title(title)\n",
" plt.colorbar().set_label('DN$^{-}$/s')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"\n",
"# 3. Create Source Catalogs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"See the [Catalog Generator Notebook](https://github.com/spacetelescope/mirage/blob/master/examples/Catalog_Generation_Tools.ipynb) for the full suite of catalog creation examples.\n",
"\n",
"There are [9 different types of source catalogs](https://mirage-data-simulator.readthedocs.io/en/latest/catalogs.html) accepted by Mirage. The types of catalogs that are accepted depend on the observing mode being simulated, as well as the type of sources in the catalogs. \n",
"\n",
"In this notebook we focus on imaging mode simulations of sidereal targets. In this case, there are three main types of catalogs that can be used:\n",
"\n",
"* Point sources\n",
"\n",
" Point source catalogs contain only point sources.\n",
"\n",
" \n",
"* Galaxies\n",
"\n",
" Galaxy catalogs contain galaxies to be added to the simulation. Galaxies are simulated as 2D Sersic profiles.\n",
"\n",
"\n",
"* \"Extended\" sources\n",
"\n",
" \"Extended\" source catalogs are a catch-all, intended for sources with more complex morphologies. In this case, the user provides a stamp image of each source, which is added to the scene.\n",
"\n",
"\n",
"Examples of these catalogs are shown below."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One important note on catalogs: every source must have its own unique index number. This includes a scenario where you are using multiple catalogs (e.g. a point source and a galaxy catalog). This is because Mirage creates a segmentation map of the scene using the index numbers. In the examples below, the 'starting_index' parameter is used to set the starting index number in a given catalog. If creating multiple catalogs, be sure to adjust indexes so that there is no overlap."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Using built-in convenience functions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Mirage contains a number of convenience functions for creating source catalogs. These functions query and create source catalogs from 2MASS, Gaia, WISE, and optionally, the Besancon Galaxy model. For each source, Mirage will use the reported magnitudes through the passbands from the various surveys, and interpolate to find the source's magnitude in the requested JWST passbands.\n",
"\n",
"Getting results from the Besancon query is a multi-step process, which is outlined in the [Catalog Generator Notebook](https://github.com/spacetelescope/mirage/blob/master/examples/Catalog_Generation_Tools.ipynb). For the purposes of this example, we will skip the Besancon query. Note that this is an optional input to the catalog creation function below.\n",
"\n",
"More description of the catalog generation conevenience tools are given in the [Catalog Creation](https://mirage-data-simulator.readthedocs.io/en/latest/catalog_creation.html) documentation page. \n",
"\n",
"The highest-level of these is the `for_proposal()` function. This function will look at the targets in your APT file and create point source and/or galaxy catalogs in the areas around those targets.\n",
"\n",
"Another useful high-level function is `get_all_catalogs()`. This function creates a point source catalog for a given central RA, Dec and box width on the sky. Under the hood, `for_proposal()` calls `get_all_catalogs()` once for each target in your APT file."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Point source catalogs created with the functions above are generated through queries to the 2MASS, Gaia, and WISE source catalogs. This will limit the sources in the resulting catalog to the magnitude limits of those three surveys. You can also optionally provide a file containing the results of a query of the Besancon Galaxy model. With this, you can add dimmer sources with a realistic source density and distribution of stellar types, but the sources will not be real stars.\n",
"\n",
"\n",
"* The catalog of extragalactic sources is created from the sources in the 3DHST catalog. As with the Besancon model results, representative galaxies will be added to the source catalog in order to create a realistic scene, but the sources will not be real. Note that in this case, the density of galaxies will match that seen in 3DHST (ie in a deep exposure)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"for_proposal()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `for_proposal()` function will generate point source and/or galaxy catalogs for the target locations contained in your APT proposal. If your proposal contains multiple targets, separate catalogs will be generated for each. The catalogs contain source magnitudes for all sources in all NIRCam/NIRISS filters specified in the proposal."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"catalog_dir = 'catalogs'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"catalog_results = create_catalog.for_proposal(xml_filename, pointing_filename,\n",
" point_source=True, extragalactic=True,\n",
" catalog_splitting_threshold=1.,\n",
" besancon_catalog_file=None,\n",
" ra_column_name='RAJ2000',\n",
" dec_column_name='DECJ2000',\n",
" out_dir=catalog_dir,\n",
" save_catalogs=True)\n",
"ptsrc_cats, gal_cats, ptsrc_filenames, gal_filenames, ptsrc_mapping, gal_mapping = catalog_results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's take a quick look at a couple of the catalogs:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"len(ptsrc_cats)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(ptsrc_filenames)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ptsrc_mapping"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"ptsrc_cats[1].table"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"ptsrc_cats[0].table"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The galaxy source catalogs are very large! We will look at the catalogs here, but not include them when creating the simulated data, in order to save time."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(gal_mapping)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"gal_cats[1].table"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"get_all_catalogs()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next let's look at a call of `get_all_catalogs()`, to show how to construct a catalog for a given area on the sky. In this case, we need to specify the JWST filters to include in the catalog. There is also an optional keyword to specify the starting index number of the catalog. This is important, because when using multiple catalogs to create a simulation, every source must have a unique index. In this case we'll assume there are no other catalogs being used, and will start with an index of 1. As with `for_proposal()`, there is an option to supply a file containing the results of a Besancon query. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"center_ra = 12.0 # degrees\n",
"center_dec = 12.0 # degrees\n",
"width = 140 # arcseconds"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sim_filters = ['F150W', 'F444W']"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cat, filt_list = create_catalog.get_all_catalogs(center_ra, center_dec, width,\n",
" besancon_catalog_file=None,\n",
" instrument='NIRCAM', filters=sim_filters,\n",
" starting_index=1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cat.table"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To use this catalog in a Mirage simulation, save it to an ascii file."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cat_filename = os.path.join(catalog_dir, 'ptsrc_from_get_all_catalogs.cat')\n",
"cat.save(cat_filename)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"galaxy_background()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Galaxies are simulated as 2D Sersic profiles. The `galaxy_background()` function uses a catalog from the [3DHST project](https://archive.stsci.edu/prepds/3d-hst/) to create catalogs of galaxies that can be simulated in Mirage. The function selects a number of galaxies from the 3DHST catalog such that the density of sources on the sky matches that in the 3DHST catalog. This is a good option for simulating a deep field observation. The `galaxy_background()` function is called by `for_proposal()`, as seen in the [example above](#for_proposal)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"center_ra = 12.0\n",
"center_dec = 12.0\n",
"v3_angle = 0. # degrees\n",
"width = 140 # arcseconds\n",
"instrument = 'nircam'\n",
"filter_list = ['F444W', 'F150W']\n",
"background_galaxy_catalog, used_seed_value = create_catalog.galaxy_background(center_ra, center_dec, v3_angle,\n",
" width, instrument, filter_list,\n",
" boxflag=False, brightlimit=14.0,\n",
" starting_index=1000)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"background_galaxy_catalog.table"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"background_galaxy_catalog.save(os.path.join(catalog_dir, 'background_galaxies_from_3DHST.cat'))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Manual creation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Point Sources"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"#### From an existing JHK catalog"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's create a small catalog containing RA, Dec, and J, H, K, and V magnitudes, plus extinction, $A_{v}$. Then we'll use a convenience function to convert it into a Mirage-formatted point source catalog."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ra_vals = [12.0, 12.001, 12.002, 12.003, 12.004]\n",
"dec_vals = [34.5, 34.5001, 34.5002, 34.5003, 34.5004]\n",
"num_stars = len(ra_vals)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"J = np.random.uniform(low=14, high=16, size=num_stars)\n",
"H = np.random.uniform(low=14, high=16, size=num_stars)\n",
"K = np.random.uniform(low=14, high=16, size=num_stars)\n",
"V = np.random.uniform(low=14, high=16, size=num_stars)\n",
"Av = np.repeat(1.0, num_stars)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"orig_cat = Table()\n",
"orig_cat['RA'] = ra_vals\n",
"orig_cat['Dec'] = dec_vals\n",
"orig_cat['J'] = J\n",
"orig_cat['H'] = H\n",
"orig_cat['K'] = K\n",
"orig_cat['V'] = V\n",
"orig_cat['Av'] = Av"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"orig_cat"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Save our \"existing\" catalog to an ascii file"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"orig_cat_file = os.path.join(catalog_dir, 'original_ptsrc_catalog.cat')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ascii.write(orig_cat, orig_cat_file, overwrite=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now convert the catalog into a Mirage-formatted catalog with magnitudes converted into the those for the filters of interest. For NIRCam, the filter needs to be specified as a filter/pupil pair. For NIRISS, only the filter name is needed."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"filters = {}\n",
"filters['nircam'] = ['F150W/CLEAR', 'F444W/CLEAR']\n",
"filters['niriss'] = ['F200W', 'F277W']"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Name of the file to save the Mirage catalog into\n",
"mirage_ptsrc_cat_file = os.path.join(catalog_dir, 'mirage_formatted_point_sources.cat')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Be sure to specify the column names in the original catalog\n",
"# that contain the RA and Dec data.\n",
"#\n",
"# Since this is the first catalog we are creating, start the index counter\n",
"# at 1. (Don't start at 0 since this will be used to create a segmentation\n",
"# map)\n",
"ptsrc_cat = create_catalog.johnson_catalog_to_mirage_catalog(orig_cat_file, filters,\n",
" ra_column_name='RA',\n",
" dec_column_name='Dec',\n",
" magnitude_system='abmag',\n",
" output_file=mirage_ptsrc_cat_file,\n",
" starting_index=1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"ptsrc_cat.table"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"#### Manually create a Mirage-formatted point source catalog from scratch"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case we create a Mirage point source catalog directly. We'll use this catalog in the simulation below, since we can tailor this catalog more easily to create a pretty picture."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Generate point source RA, Dec values that cover the field of view\n",
"# for all detectors. \n",
"min_ra = 11.980270819703372\n",
"max_ra = 12.050540819703372\n",
"min_dec = 11.965394574641623\n",
"max_dec = 12.035664574641622\n",
"delta_ra = max_ra - min_ra\n",
"delta_dec = max_dec - min_dec"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Generate a list of RA, Dec pairs\n",
"num_stars = 450\n",
"random_number_generator = np.random.RandomState(2021)\n",
"ra = random_number_generator.rand(num_stars) * delta_ra + min_ra\n",
"dec = random_number_generator.rand(num_stars) * delta_dec + min_dec\n",
"\n",
"# Create a list of magnitudes for two filters. Let's keep the magnitudes\n",
"# between 14 and 20, just to keep things easily visible in this short\n",
"# exposure\n",
"mag_rand_num_gen = np.random.RandomState(1066)\n",
"mags1 = mag_rand_num_gen.rand(num_stars) * 6 + 14.\n",
"mags2 = mag_rand_num_gen.rand(num_stars) * 6 + 14."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create a PointSourceCatalog object and supply the RA and Dec values.\n",
"# We won't use this catalog for our simulations. \n",
"ptsrc = catalog_generator.PointSourceCatalog(ra=ra, dec=dec, starting_index=1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Now add magnitude columns for each filter\n",
"ptsrc.add_magnitude_column(mags1, instrument='nircam',\n",
" filter_name='f150w', magnitude_system='abmag')\n",
"ptsrc.add_magnitude_column(mags2, instrument='nircam',\n",
" filter_name='f444w', magnitude_system='abmag')\n",
"ptsrc.add_magnitude_column(mags2, instrument='niriss',\n",
" filter_name='f200w', magnitude_system='abmag')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"manually_generated_ptsrc_catalog = os.path.join(catalog_dir, 'manually_generated_ptsrc.cat')\n",
"ptsrc.save(manually_generated_ptsrc_catalog)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"ptsrc.table"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"See the [Catalog Generator Notebook](https://github.com/spacetelescope/mirage/blob/master/examples/Catalog_Generation_Tools.ipynb) for more examples of creating source catalogs using queries to 2MASS/GAIA/WISE/Besancon."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Galaxy catalog from scratch"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case we create a Mirage galaxy source catalog directly. We'll use this catalog in the simulation below, since the galaxy catalog created above is so large, and would take longer to run."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"num_galaxies = 20\n",
"gal_rand_num_gen = np.random.RandomState(1564)\n",
"ra_galaxy_vals = gal_rand_num_gen.rand(num_galaxies) * delta_ra + min_ra\n",
"dec_galaxy_vals = gal_rand_num_gen.rand(num_galaxies) * delta_dec + min_dec"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"radius = np.random.uniform(low=0.06, high=0.5, size=num_galaxies) # arcsec\n",
"ellip = np.random.uniform(low=0., high=0.8, size=num_galaxies)\n",
"posang = np.random.uniform(low=0, high=359, size=num_galaxies) # degrees\n",
"sersic = np.random.uniform(low=1.0, high=4.0, size=num_galaxies)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Manually add a galaxy at a known location so we can examine it later\n",
"ra_galaxy_vals[-1] = 12.007490819703373\n",
"dec_galaxy_vals[-1] = 11.992614574641623"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Tweak the properties of the galaxy in the known location\n",
"radius[-1] = 0.1\n",
"ellip[-1] = 0.7\n",
"posang[-1] = 35.\n",
"sersic[-1] = 1.5"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Since we already have the point source catalog, we know the minimum index\n",
"# value that we can use in order to ensure that all sources have a unique index.\n",
"# Set the minimum index number to be one larger than the number of point sources.\n",
"gal_cat = catalog_generator.GalaxyCatalog(ra=ra_galaxy_vals, dec=dec_galaxy_vals,\n",
" ellipticity=ellip,\n",
" radius=radius,\n",
" sersic_index=sersic,\n",
" position_angle=posang,\n",
" radius_units='arcsec',\n",
" starting_index=num_stars+1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"gal_mag_f150w = mag_rand_num_gen.rand(num_galaxies) + 15.\n",
"gal_mag_f444w = mag_rand_num_gen.rand(num_galaxies) + 15.5"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"gal_cat.add_magnitude_column(gal_mag_f150w, instrument='nircam', filter_name='f150w',\n",
" magnitude_system='abmag')\n",
"gal_cat.add_magnitude_column(gal_mag_f444w, instrument='nircam', filter_name='f444w',\n",
" magnitude_system='abmag')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"manually_generated_galaxy_catalog = os.path.join(catalog_dir, 'manually_generated_galaxies.cat')\n",
"gal_cat.save(manually_generated_galaxy_catalog)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"gal_cat.table"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### \"Extended\" source catalog from scratch"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\"Extended\" sources can be used for sources where you have a fits file containing a stamp image of your source. In this way, Mirage can simulate sources with more complex morphologies than simple point sources and Sersic profiles. Here, we'll create an extended source catalog with one source, in order to show its use. We'll set it up so the extended source falls onto the B4 detector."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Download a stamp image to use for this example\n",
"box_file = 'https://stsci.box.com/shared/static/lkfn6oz03pbyorka0644x2x8344q3oyg.fits'\n",
"stamp_file = 'extended_stamp.fits'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"urllib.request.urlretrieve(box_file, stamp_file)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's take a quick look at this stamp image. This will be the object we are adding to the scene via the extended source catalog."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"obj = fits.getdata(stamp_file)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show(obj, 'Extended source stamp image', min=0, max=10)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"extended_ra = [12.0088178]\n",
"extended_dec = [11.9911822]\n",
"extended_stamp_file = ['extended_stamp.fits']\n",
"extended_pa = [0.]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Note that the highest index number in the galaxy catalog is 60, so \n",
"# set the starting_index to something higher.\n",
"extended_cat = catalog_generator.ExtendedCatalog(filenames=extended_stamp_file,\n",
" ra=extended_ra,\n",
" dec=extended_dec,\n",
" position_angle=extended_pa,\n",
" starting_index=99999)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For extended sources, you can either specify a magnitude, as with the point sources/galaxies, or you can specify the magnitude as 'None'. In the latter case, the stamp image is interpreted as being in units of counts per second, and is added directly to the simulation without any scaling."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"extended_mag_f150w = [14.]\n",
"extended_mag_f444w = [16.]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"extended_cat.add_magnitude_column(extended_mag_f444w, instrument='nircam', \n",
" filter_name='f444w', magnitude_system='abmag')\n",
"extended_cat.add_magnitude_column(extended_mag_f150w, instrument='nircam', \n",
" filter_name='f150w', magnitude_system='abmag')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"extended_cat.table"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"manually_generated_extended_catalog = os.path.join(catalog_dir, 'manually_generated_extended.cat')\n",
"extended_catalog_file = os.path.join(catalog_dir, 'manually_generated_extended.cat')\n",
"extended_cat.save(manually_generated_extended_catalog)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Back to the [Table of contents](#toc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"\n",
"# Generating input yaml files"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For convenience, observing programs with multiple pointings \n",
"and detectors can be simulated starting with the program's \n",
"APT file. The xml and pointings files must be exported from \n",
"APT, and are then used as input to the *yaml_generator*, which will\n",
"generate a series of yaml input files. The [yaml generator documentation page](https://mirage-data-simulator.readthedocs.io/en/latest/yaml_generator.html) explains the\n",
"creation of yaml files in more detail."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Optional user inputs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"See Mirage's [Mirage's yaml_generator documentation](https://mirage-data-simulator.readthedocs.io/en/latest/yaml_generator.html#additional-yaml-generator-inputs \"Yaml Generator Inputs\")\n",
"for details on the formatting options for the inputs listed below. The formats will vary based on the complexity of your inputs."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Catalogs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Catalogs are organized by target name. If starting from an APT file, these must be the target names in your proposal. Mirage will then map the specified catalogs to the appropriate observations that use each target.\n",
"\n",
"There are several ways to specify the catalogs when calling the `yaml_generator`. These are detailed on the [yaml generation documentation page](https://mirage-data-simulator.readthedocs.io/en/latest/yaml_generator.html#source-catalogs). In this notebook, we will specify one catalog of each type for each target in the APT file."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```python\n",
"# When using the for_proposal() or get_all_catalogs() convenience functions, you can\n",
"# populate catalog names directly from the output names. \n",
"\n",
"cat_dict = {'STARFIELD': {'point_source': os.path.join(catalog_dir, ptsrc_mapping['001']),\n",
" 'extended': manually_generated_extended_catalog\n",
" },\n",
" 'EXTRAGALACTIC': {'point_source': os.path.join(catalog_dir, ptsrc_mapping['002']),\n",
" 'galaxy': os.path.join(catalog_dir, gal_mapping['002'])\n",
" }\n",
" }\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case, let's use our manually-created catalogs for the STARFIELD target, where we were able to control the number of sources in the field of view. Our manually-created catalogs do not extend to the EXTRAGALACTIC source, so for those we continue to use the outputs from for_proposal()."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cat_dict = {'STARFIELD': {'point_source': manually_generated_ptsrc_catalog,\n",
" 'galaxy': manually_generated_galaxy_catalog, \n",
" 'extended': manually_generated_extended_catalog\n",
" },\n",
" 'EXTRAGALACTIC': {'point_source': os.path.join(catalog_dir, ptsrc_mapping['002']),\n",
" 'galaxy': os.path.join(catalog_dir, gal_mapping['002'])\n",
" }\n",
" }"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Pipeline reference files"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Set reference file values. Setting this to 'crds_full_name' when calling the `yaml_generator` will cause the yaml_generator to search for and download needed calibration reference files (commonly referred to as CRDS reference files) when the yaml_generator is run. The names of these CRDS reference files will then be placed in the appropriate entries of the yaml files. This option can be useful if you want to be able to guarantee the use of the same reference files no matter when the yaml file is used to create a simulation. \n",
" \n",
"Setting this to 'crds' will put placeholders (the string \"crds\") in the yaml files' entries for CRDS reference files. In this case, the CRDS reference files are not identified and downloaded until the simulated data are created from the yaml file. With this method, you are guaranteed to use the latest CRDS reference files when creating the simulated data, even if your yaml files are old."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"reffile_defaults = 'crds'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Cosmic rays"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can control the library and cosmic ray rate if desired. If you omit this line, Mirage will use the default values. More details on rates are given on the [Observation Generation documentation page](https://mirage-data-simulator.readthedocs.io/en/latest/observation_generator.html#add-cosmic-rays) as well as the [yaml generation documentation page](https://mirage-data-simulator.readthedocs.io/en/latest/yaml_generator.html#cosmic-ray-rates)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cosmic_rays = {'library': 'SUNMAX', 'scale': 1.0}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Background"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This option controls the background signal in the simulations. Mirage uses the same definitions as ETC for the options: 'low', 'medium', and 'high'. Or, if you set background to a number, Mirage will assume that this is the background in counts per second per pixel.\n",
"\n",
"Another option for the background is to specify that you want the background level associated with a particular date. If the `dateobs_for_background` parameter is set to True in the call to the yaml_generator, then any background value given here is ignored and the background level will be calculated using the `jwst_backgrounds` tool.\n",
"\n",
"More details are given on the [yaml generation documentation page](https://mirage-data-simulator.readthedocs.io/en/latest/yaml_generator.html#background-specification)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Set the background for all observations\n",
"background = 'medium'\n",
"\n",
"# Give a different background value for each observation,\n",
"# where the keys here are observation numbers from the APT file\n",
"#background = {'001': 'high', '002': 'medium'}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Telescope Roll Angle"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can set the telescope roll angle on a per proposal or per observation basis. Note that this is the roll angle about the V1 axis in degrees east of North. If you omit this parameter, the default is a roll angle of 0. From this value, Mirage will calculate the local roll angle of the detector to be simulated.\n",
"\n",
"More details are shown on the [yaml generation documentation page](https://mirage-data-simulator.readthedocs.io/en/latest/yaml_generator.html#roll-angle)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Set one roll angle for all APT observations:\n",
"pav3 = 0.\n",
"\n",
"# Or set a different roll angle for each APT observation. In this way you\n",
"# can simulate different epochs.\n",
"#roll_angle = {'001': 34.5, '002': 154.5}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Date"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Set the observation date to use for the data. A single date can be given for the proposal, or separate dates can be provided for each observation within the proposal. This information is placed in the headers of the output files. If `dateobs_for_background` is set to True in the call to the `yaml_generator`, then the date will be used to calculate the background signal.\n",
"\n",
"See the [yaml generator documentation page](https://mirage-data-simulator.readthedocs.io/en/latest/yaml_generator.html#observation-dates) for more details."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Set one date for all APT observations\n",
"dates = '2022-10-31'\n",
"\n",
"# Specify a different date for each APT observation\n",
"#dates = {'001': '2022-06-25', '002': '2022-11-15'}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Ghosts"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For NIRISS simulations, users can add optical ghosts to the data. By default, ghosts will be added for point sources only. Ghosts can also be added for galaxy or extended targets if you have a stamp image for each source. See the [documentation for adding ghosts](https://mirage-data-simulator.readthedocs.io/en/latest/ghosts.html)\n",
"for details. For NIRCam simulations, such as those created in this notebook, the addition of ghosts is not supported, and Mirage will ignore the keywords below, but we include them here for completeness."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ghosts = False\n",
"convolve_ghosts = False"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Run the yaml_generator"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This will create a collection of yaml files that will be used as inputs when creating the simulated data. There will be one yaml file for each detector and exposure, so there can be quite a few files created if your program has lots of exposures or dithers. A more complete description of a call to the `yaml_generator` is given on the [yaml generation documentation page](https://mirage-data-simulator.readthedocs.io/en/latest/yaml_generator.html#run-the-yaml-generator)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Set the directory into which the yaml files will be written\n",
"output_dir = './yaml_files/'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# You can also set a separate directory where the simulated data\n",
"# will eventually be saved to\n",
"simulation_dir = './sim_data/'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can specify the data reduction state of the Mirage outputs.\n",
"Options are 'raw', 'linear', or 'linear, raw'. \n",
"\n",
"If 'raw' is specified, the output is a completely uncalibrated file, with a filename ending in \"uncal.fits\"\n",
"\n",
"If 'linear' is specified, the output is a file with linearized signals, ending in \"linear.fits\". This is equivalent to having been run through the dq_init, saturation flagging, superbias subtraction, reference pixel subtraction, and non-linearity correction steps of the calibration pipeline. Note that this product does not include dark current subtraction.\n",
"\n",
"If 'linear, raw', both outputs are saved.\n",
"\n",
"In order to fully process the Mirage output with the default steps used by the pipeline, it is **recommended to use the 'raw' output and run the entire calibration pipeline after Mirage has created all the data**."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"datatype = 'raw, linear'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"# Run the yaml generator\n",
"yam = yaml_generator.SimInput(input_xml=xml_filename, pointing_file=pointing_filename,\n",
" catalogs=cat_dict, cosmic_rays=cosmic_rays,\n",
" background=background, roll_angle=pav3,\n",
" dates=dates, reffile_defaults=reffile_defaults,\n",
" add_ghosts=ghosts, convolve_ghosts_with_psf=convolve_ghosts,\n",
" verbose=True, output_dir=output_dir,\n",
" simdata_output_dir=simulation_dir,\n",
" datatype=datatype)\n",
"yam.create_inputs()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"yfiles = sorted(glob(os.path.join(output_dir,'jw*.yaml')))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"yfiles"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Back to the [Table of contents](#toc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Examine a yaml input file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The yaml input file contains all of the parameters and values that Mirage needs in order to simulate one exposure from one detector. Keep in mind that the yaml generator above is a convenience function for quickly generating the yaml files associated with a particular proposal. If desired, you can always create your own yaml files, or take an existing yaml file and tweak it in order to customize your simulation.\n",
"\n",
"The Mirage documentation also provides an [example of a yaml file](https://mirage-data-simulator.readthedocs.io/en/latest/example_yaml.html), complete with descriptions for all\n",
"paramteres.\n",
"\n",
"Entries listed as 'config' have default files that are present in the \n",
"config directory of the repository. The scripts are set up to \n",
"automatically find and use these files. The user can replace 'config'\n",
"with a filename if they wish to override the default.\n",
"\n",
"In general, if 'None' is placed in a field, then the step that uses\n",
"that particular file will be skipped.\n",
"\n",
"Note that the linearized_darkfile entry overrides the dark entry, unless\n",
"linearized_darkfile is set to None, in which case the dark entry will be\n",
"used.\n",
"\n",
"Use of a valid readout pattern in the readpatt entry will cause the \n",
"simulator to look up the values of nframe and nskip and ignore the \n",
"values given in the yaml file."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's take a quick look at one of the yaml files that were created above."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Choose one of the yaml files just created\n",
"yamlfile = './yaml_files/jw98765001001_01101_00001_nrcb4.yaml'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"with open(yamlfile) as f:\n",
" yaml_data = yaml.load(f, Loader=yaml.FullLoader)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"yaml_data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Back to the [Table of contents](#toc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"\n",
"# Create simulated data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Under the hood, the `Mirage` simulator is broken up into three basic stages:\n",
"\n",
"1. **Creation of a \"seed image\".**
\n",
" This is generally a noiseless countrate image that contains signal\n",
" only from the astronomical sources to be simulated. Currently, the \n",
" mirage package contains code to produce a seed image starting\n",
" from object catalogs.
\n",
" \n",
"2. **Dark current preparation.**
\n",
" The simualted data will be created by adding the simulated sources\n",
" in the seed image to a real dark current exposure. This step\n",
" converts the dark current exposure to the requested readout pattern\n",
" and subarray size requested by the user.
\n",
" \n",
"3. **Observation generation.**
\n",
" This step converts the seed image into an exposure of the requested\n",
" readout pattern and subarray size. It also adds cosmic rays and \n",
" Poisson noise, as well as other detector effects (IPC, crosstalk, etc).\n",
" This exposure is then added to the dark current exposure from step 2.
\n",
" \n",
"For imaging mode observations, these steps are wrapped by the `imaging_simulator.py` module, as shown below."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Call the imaging simulator"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The imaging_simulator.ImgSim class is a wrapper around the three main steps of the simulator (detailed in the [Running simulator steps independently](#run_steps_independently) section below). This convenience function is useful when creating simulated imaging mode data. WFSS data will need to be run in a slightly different way. See the WFSS example notebook for details."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"# Run all steps of the imaging simulator for yaml file #1\n",
"img_sim = imaging_simulator.ImgSim()\n",
"img_sim.paramfile = yamlfile\n",
"img_sim.create()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Examine the Output"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Noiseless Seed Image"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This image is an intermediate product. It contains only the signal from the astronomical sources and background. There are no detector effects, nor cosmic rays added to this count rate image."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# First, look at the noiseless seed image\n",
"show(img_sim.seedimage,'Seed Image', max=50)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"# See the galaxy source\n",
"show(img_sim.seedimage[975:1075, 975:1075],'Galaxy in Seed Image', max=2000)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# See the extended source\n",
"show(img_sim.seedimage[800:900, 820:920],'Extended Source in Seed Image', max=2000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Raw (uncal) file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is the \"final\" output for Mirage. The uncal file contains raw, uncalibrated data, and is in a format that matches the level 1b data that will be produced by JWST. This file may be run through the JWST calibration pipeline just as if it were real JWST data. Note that Mirage's creation of the raw and linearized output exposures is controlled through the `datatype` parameter in the yaml generator."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Examine the raw output. First a single group, which is dominated by noise and detector artifacts. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"raw_basename = os.path.basename(yamlfile).replace('.yaml', '_uncal.fits')\n",
"raw_file = os.path.join(simulation_dir, raw_basename)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"hdulist = fits.open(raw_file)\n",
"hdulist.info()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"raw_data = hdulist['SCI'].data\n",
"raw_header = hdulist[0].header\n",
"hdulist.close()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(raw_data.shape)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show(raw_data[0, -1, :, :], \"Final Group\", max=15000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Many of the instrumental artifacts can be removed by looking at the difference between two groups. Raw data values are integers, so first make the data floats before doing the subtraction. Here we will look at the difference between the last and the first groups."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show(1. * raw_data[0, -1, :, :] - 1. * raw_data[0, 0, :, :],\n",
" \"Last Minus First Group\", max=2000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This raw data file is now ready to be run through the [JWST calibration pipeline](https://jwst-pipeline.readthedocs.io/en/stable/) from the beginning. If dark current subtraction is not important for you, you can use Mirage's linear output, skip some of the initial steps of the pipeline, and begin by running the [Jump detection](https://jwst-pipeline.readthedocs.io/en/stable/jwst/jump/index.html?highlight=jump) and [ramp fitting](https://jwst-pipeline.readthedocs.io/en/stable/jwst/ramp_fitting/index.html) steps."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Linearized exposure"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Another optional output is a version of the raw exposure above that contains linearized data. In this case, the data are saved in a state equvalent to that output from the linearization step of the calibration pipeline. Visual examination of linearized data is often easier than that of raw data, because the linearized data has had bias drifts removed through the use of reference pixels. However, the signal to noise in this file will be lower than that in the seed image, since the linearized exposure does contain noise as well as containing full integrations with multiple groups, rather than the line-fit slope image present in the seed image."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"linearized_file = os.path.join(simulation_dir, 'jw98765001001_01101_00001_nrcb4_linear.fits')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"linearized_data = fits.getdata(linearized_file)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's look at the signal difference between the raw file and the linearized file for a pixel within a source."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"linearized_data[0, :, 1025, 1025]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"raw_data[0, :, 1025, 1025]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"times = raw_header['TGROUP'] * np.arange(raw_header['NGROUPS'])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"slope = (raw_data[0, 1, 1025, 1025] - raw_data[0, 0, 1025, 1025]) / (times[1] - times[0])\n",
"straight_line = slope * times"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"min_lin = linearized_data[0, 0, 1025, 1025]\n",
"min_raw = raw_data[0, 0, 1025, 1025]\n",
"\n",
"f, a = plt.subplots(figsize=(8,8))\n",
"a.plot(times, linearized_data[0, :, 1025, 1025] - min_lin, 'o-', color='red', label='Linearized')\n",
"a.plot(times, raw_data[0, :, 1025, 1025] - min_raw, 'o-', color='blue', label='Raw')\n",
"a.plot(times, straight_line, color='black', linestyle=(0, (5, 10)), label='Extrapolated from raw')\n",
"a.set_xlabel('Time (sec)')\n",
"a.set_ylabel('Signal (ADU)')\n",
"a.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Back to the [Table of contents](#toc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"\n",
"# Running simulation steps independently"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The steps detailed in this section are wrapped by the `imaging_simulator` mentioned above. General users will not need to worry about the details of these three steps."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## First generate the \"seed image\" \n",
"\n",
"This is generally a 2D noiseless countrate image that contains only simulated astronomical sources.\n",
"\n",
"A seed image is generated based on a `.yaml` file that contains all the necessary parameters for simulating data. For this exercise, use the same yaml file that was used in the [Create Simulated Data](#run_steps_together) section as input."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"cat = catalog_seed_image.Catalog_seed()\n",
"cat.paramfile = yamlfile\n",
"cat.make_seed()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Look at the seed image"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show(cat.seedimage,'Seed Image',max=50)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"While the seed image makes for a pretty picture, and is useful as a sanity check that your objects contain the correct signals and are at the correct locations, it is far from a complete simulation. It contains no noise (Poisson, readnoise, 1/f noise, etc). Also, there is no WCS information attached to the seed image file, so it cannot be run through the calibration pipeline."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"seed_image_hdulist = fits.open(cat.seed_file)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"seed_image_hdulist.info()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that the seed image file contains an extension with a segmentation map. For imaging mode simulations, Mirage does not use this output product at all. For WFSS simulations, the segmentation map controls which pixels are dispersed when creating the seed image. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show(seed_image_hdulist[2].data,'Segmentation Map', max=400)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"seed_image_header = seed_image_hdulist[0].header"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"seed_image_header"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Other output products"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Catalog of sources present on the detector"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"img_point_source_cat = os.path.join(simulation_dir, 'jw98765001001_01101_00001_nrcb4_uncal_pointsources.list')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"img_point_sources = ascii.read(img_point_source_cat, format='commented_header', header_start=2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"img_point_sources"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"img_gal_source_cat = os.path.join(simulation_dir, 'jw98765001001_01101_00001_nrcb4_uncal_galaxySources.list')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"img_gal_sources = ascii.read(img_gal_source_cat)\n",
"img_gal_sources"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"img_ext_source_cat = os.path.join(simulation_dir, 'jw98765001001_01101_00001_nrcb4_uncal_extendedsources.list')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"img_ext_sources = ascii.read(img_ext_source_cat, format='commented_header', header_start=2)\n",
"img_ext_sources"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Catalog of cosmic rays added to the exposure"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cr_file = os.path.join(simulation_dir, 'jw98765001001_01101_00001_nrcb4_uncal_cosmicrays.list')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cr_table = ascii.read(cr_file)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"cr_table"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Log files"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Log files for completed Mirage runs will be saved in the `mirage_logs` subdirectory under the working directory."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"logfiles = sorted(glob('mirage_logs/*.log'))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"with open(logfiles[-1]) as obj:\n",
" log = obj.readlines()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"log"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The log file for the latest Mirage run will also be stored in the current working directly as `mirage_latest.log`. If Mirage crashes, this is the log file you should examine, as Mirage copies `mirage_latest.log` into the `mirage_logs` directory at the completion of the run."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"with open('mirage_latest.log') as obj:\n",
" latest_log = obj.readlines()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"latest_log"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Prepare the dark current exposure\n",
"This will serve as the base of the simulated data.\n",
"This step will linearize the dark current (if it \n",
"is not already), and reorganize it into the \n",
"requested readout pattern and number of groups."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"d = dark_prep.DarkPrep()\n",
"d.paramfile = yamlfile\n",
"d.prepare()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Look at the dark current \n",
"For this, we will look at an image of the final group\n",
"minus the first group"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"exptime = d.linDark.header['NGROUPS'] * cat.frametime\n",
"diff = (d.linDark.data[0,-1,:,:] - d.linDark.data[0,0,:,:]) / exptime\n",
"show(diff,'Dark Current Countrate',max=0.001)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"darkfile = 'sim_data/jw98765001001_01101_00001_nrcb4_uncal_linear_dark_prep_object.fits'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dark_header = fits.getheader(darkfile)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dark_header"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Create the final exposure\n",
"Turn the seed image into a exposure of the \n",
"proper readout pattern, and combine it with the\n",
"dark current exposure. Cosmic rays and other detector\n",
"effects are added. \n",
"\n",
"The output can be either this linearized exposure, or\n",
"a 'raw' exposure where the linearized exposure is \n",
"\"unlinearized\" and the superbias and \n",
"reference pixel signals are added, or the user can \n",
"request both outputs. This is controlled from\n",
"within the yaml parameter file."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"obs = obs_generator.Observation()\n",
"obs.linDark = d.prepDark\n",
"obs.seed = cat.seedimage\n",
"obs.segmap = cat.seed_segmap\n",
"obs.seedheader = cat.seedinfo\n",
"obs.paramfile = yamlfile\n",
"obs.create()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Examine the final output image\n",
"Again, we will look at the last group minus the first group"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"obs.linear_output"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"with fits.open(obs.linear_output) as h:\n",
" lindata = h[1].data\n",
" header = h[0].header"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"exptime = header['EFFINTTM']\n",
"diffdata = (lindata[0,-1,:,:] - lindata[0,0,:,:]) / exptime\n",
"show(diffdata,'Simulated Data',min=0,max=50)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Show on a log scale, to bring out the presence of the dark current.\n",
"# Noise in the CDS image makes for a lot of pixels with values < 0,\n",
"# which makes this kind of an ugly image. Add an offset so that\n",
"# everything is positive and the noise is visible\n",
"offset = 2.\n",
"plt.figure(figsize=(12,12))\n",
"plt.imshow(np.log10(diffdata+offset),clim=(0.001,np.log10(50)), origin='lower')\n",
"plt.title('Simulated Data')\n",
"plt.colorbar().set_label('DN$^{-}$/s')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"\n",
"## Simulating Multiple Exposures"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Each yaml file will simulate an exposure for a single pointing using a single detector. Here, let's simulate the data from all 5 B module detectors for a single pointing. To save time, we'll create only the stamp images."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"first_pointing_yamls = ['./yaml_files/jw98765001001_01101_00001_nrcb1.yaml',\n",
" './yaml_files/jw98765001001_01101_00001_nrcb2.yaml',\n",
" './yaml_files/jw98765001001_01101_00001_nrcb3.yaml',\n",
" './yaml_files/jw98765001001_01101_00001_nrcb4.yaml',\n",
" './yaml_files/jw98765001001_01101_00001_nrcb5.yaml']"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"# This cell will take a minute or two to run\n",
"for yfile in first_pointing_yamls:\n",
" cat = catalog_seed_image.Catalog_seed()\n",
" cat.paramfile = yfile\n",
" cat.make_seed()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"-------\n",
"If you want to skip running the cell above, but still want to display the results, you can download the seed images in the cells below."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"box_links = ['https://stsci.box.com/shared/static/r3kw8jr3l7swllspu4myqg9f4eqmvwjy.fits',\n",
" 'https://stsci.box.com/shared/static/evoq9pyu709ssoljt6my3qvu4u7eagwi.fits',\n",
" 'https://stsci.box.com/shared/static/t4ma9r2y4thzq4jstk9z99pfu0cpgwy2.fits',\n",
" 'https://stsci.box.com/shared/static/m4t472dj9uxgtna1jd551o7a86laurwt.fits',\n",
" 'https://stsci.box.com/shared/static/hkpegevuggmz42zhj2x89s6cu93xihto.fits']\n",
"deepfield_filenames = ['jw98765001001_01101_00001_nrcb1_uncal_F150W_CLEAR_final_seed_image.fits',\n",
" 'jw98765001001_01101_00001_nrcb2_uncal_F150W_CLEAR_final_seed_image.fits',\n",
" 'jw98765001001_01101_00001_nrcb3_uncal_F150W_CLEAR_final_seed_image.fits',\n",
" 'jw98765001001_01101_00001_nrcb4_uncal_F150W_CLEAR_final_seed_image.fits',\n",
" 'jw98765001001_01101_00001_nrcb5_uncal_F444W_CLEAR_final_seed_image.fits']"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for link, filename in zip(box_links, deepfield_filenames):\n",
" urllib.request.urlretrieve(link, filename)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"------\n",
"Let's arrange the seed images following the NIRCam detector layout, and view them all together. Note that seed images do not contain WCS information, so they cannot be arranged by WCS (e.g. in ds9)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"seed_img_files = []\n",
"for yfile in first_pointing_yamls:\n",
" changedir = yfile.replace('./yaml_files', './sim_data')\n",
" detector = changedir.split('_')[-1]\n",
" if 'b5' in detector:\n",
" filtername = 'F444W'\n",
" else:\n",
" filtername = 'F150W'\n",
" seedfile = changedir.replace('.yaml', '_uncal_{}_CLEAR_final_seed_image.fits'.format(filtername))\n",
" seed_img_files.append(seedfile)\n",
"seed_img_files"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Read in all of the seed images\n",
"b1 = fits.getdata(seed_img_files[0])\n",
"b2 = fits.getdata(seed_img_files[1])\n",
"b3 = fits.getdata(seed_img_files[2])\n",
"b4 = fits.getdata(seed_img_files[3])\n",
"b5 = fits.getdata(seed_img_files[4])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def display_sw(one, two, three, four, min=0, max=1000):\n",
" fig, ax = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(10, 10))\n",
" ax[0, 0].imshow(three, clim=(min, max), origin='lower')\n",
" ax[0, 1].imshow(one, clim=(min, max), origin='lower')\n",
" ax[1, 0].imshow(four, clim=(min, max), origin='lower')\n",
" ax[1, 1].imshow(two, clim=(min, max), origin='lower')\n",
" # Hide x labels and tick labels for top plots and y ticks for right plots.\n",
" for a in ax.flat:\n",
" a.label_outer()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"display_sw(b1, b2, b3, b4, min=0, max=1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show(b5, 'LW Channel', max=10)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Back to the [Table of contents](#toc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"-----\n",
"\n",
"# Simulate deep field exposure"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this section, we create a simluation of a deep-field-like exposure, using the background galaxy catalog produced from the `background_galaxy()` function. Creating and adding galaxies to the seed image takes longer than adding point sources, so this cell will take longer than the simulations above."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"deep_field_yaml = './yaml_files/jw98765002001_01101_00001_nrcb5.yaml'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# NOTE: this cell may take a few minutes to run\n",
"# If you don't want to wait for it, you can download the\n",
"# resulting seed image using the cell below.\n",
"\n",
"deepfield_sim = imaging_simulator.ImgSim()\n",
"deepfield_sim.paramfile = deep_field_yaml\n",
"deepfield_sim.create()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"-----\n",
"If you skipped running the imaging simulator above, you can download the seed image for the deep field observation here"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"box_link = 'https://stsci.box.com/shared/static/6z1yu5pju1fcbfld4u0e6nf34io2sr7g.fits'\n",
"deepfield_file = 'jw98765002001_01101_00001_nrcb5_uncal_F444W_CLEAR_final_seed_image.fits'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"urllib.request.urlretrieve(box_link, deepfield_file)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"deepfield_sim = catalog_seed_image.Catalog_seed()\n",
"deepfield_sim.seedimage = fits.getdata(deepfield_file)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"-----"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show(deepfield_sim.seedimage,'Seed Image', max=2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show(deepfield_sim.seedimage[0:500, 250:750],'Seed Image', max=2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Back to the [Table of contents](#toc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"# 10. Calibrate the data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The \"raw\" outputs from Mirage are equivalent to Level-1b data that will be returned from JWST. These files contain the suffix \"_uncal.fits\". \n",
"\n",
"You can now proceed with calibration using the JWST data calibration pipeline, just as with real data. For imaging mode data such as those produced here, a previous [JWEbbinar contains notebooks](https://github.com/spacetelescope/jwebbinar_prep/tree/main/imaging_mode) showing how to run the calibration pipeline.\n",
"\n",
"[Pipeline documentation](https://jwst-pipeline.readthedocs.io/en/latest/) is also available through readthedocs."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Back to the [Table of contents](#toc)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "JWebbinar",
"language": "python",
"name": "jwebbinar"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}