diff --git a/notebooks/stips/requirements.txt b/notebooks/stips/requirements.txt index 6f22872..37e37d1 100644 --- a/notebooks/stips/requirements.txt +++ b/notebooks/stips/requirements.txt @@ -1,2 +1,4 @@ stips -astropy \ No newline at end of file +astropy +matplotlib +esutil # --no-cache-dir option recommended on pip install \ No newline at end of file diff --git a/notebooks/stips/stips.ipynb b/notebooks/stips/stips.ipynb index 5df72d3..2d65275 100644 --- a/notebooks/stips/stips.ipynb +++ b/notebooks/stips/stips.ipynb @@ -25,6 +25,7 @@ { "cell_type": "markdown", "metadata": { + "jp-MarkdownHeadingCollapsed": true, "slideshow": { "slide_type": "slide" } @@ -36,7 +37,7 @@ "\n", "Though it trades some accuracy in order to capture the full array of detectors simulations – see the Pandeia Tutorial for higher accuracy simulations of smaller areas – STIPS simulations do take Roman's exposure-level pipeline (\"Level 2\") into account. This means generated scenes come with readouts of many calibration residuals. Scenes are also returned with Poisson and readout noise estimates and can incorporate instrumental distortion.\n", "\n", - "This notebook is a starter guide to simulating and manipulating scenes with STIPS. STIPS requires separate reference data both for itself and for some of its dependencies, so be sure you've followed the complete installation instructions before attempting to run the notebook." + "This notebook is a starter guide to simulating and manipulating scenes with STIPS. STIPS requires separate reference data both for itself and for some of its dependencies, so if you're running this notebook locally, be sure you've followed [the complete installation instructions](https://stips.readthedocs.io/en/latest/installation.html) before attempting to run it." ] }, { @@ -45,7 +46,7 @@ "source": [ "## Imports\n", "\n", - "Besides the STIPS-related imports, `astropy.io.fits` will help write a FITS table on the fly." + "Besides the STIPS-related imports, the `matplotlib` imports will help visualize simulated images and `astropy.io.fits` will help write a FITS table on the fly." ] }, { @@ -54,14 +55,13 @@ "metadata": {}, "outputs": [], "source": [ + "import matplotlib as mpl\n", + "import matplotlib.pyplot as plt\n", "import stips\n", "\n", "from astropy.io import fits\n", "from stips.observation_module import ObservationModule\n", - "from stips.scene_module import SceneModule\n", - "\n", - "# USE PYTHON 3.11\n", - "# GET EXAMPLE FILE IN LAST TUTORIAL (PSFs + adding sources)" + "from stips.scene_module import SceneModule" ] }, { @@ -70,9 +70,7 @@ "source": [ "### Environment report\n", "\n", - "To verify the existing STIPS installation alongside its associated data files and dependencies, run the cell below.\n", - "\n", - "[SAY SHOULD MATCH https://stips.readthedocs.io/en/latest/installation.html#stips-requirements ??]" + "To verify the existing STIPS installation alongside its associated data files and dependencies, run the cell below. (Find the current software requirements [in the STIPS documentation](https://stips.readthedocs.io/en/latest/installation.html#stips-requirements).)" ] }, { @@ -102,38 +100,31 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "* Installation test (overview) -- included in Intro/Import section\n", - "* Basic STIPS usage (tutorial)\n", - "* Generating scene from user-based catalog (tutorial)\n", - "* Modifying observation parameters (tutorial)\n", - "* PSFs and adding sources" + "### Basic STIPS usage\n", + "\n", + "At its most fundamental level, STIPS takes a dictionary of observation and instrument parameters and a FITS source catalog in order to return a simulated image.\n", + "\n", + "In this example, we build a dictionary specifying that the image is taken with the F129 filter, comes from the WFI01 detector, and has an exposure time of 300 seconds." ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "### Basic STIPS usage" + "obs = {'instrument': 'WFI', 'filters': ['F129'], 'detectors': 1,\n", + " 'background': 'pandeia', 'observations_id': 42, 'exptime': 300,\n", + " 'offsets': [{'offset_id': 1, 'offset_centre': False, 'offset_ra': 0.0,\n", + " 'offset_dec': 0.0, 'offset_pa': 0.0}]\n", + " }" ] }, { - "cell_type": "raw", + "cell_type": "markdown", "metadata": {}, "source": [ - "# Build observation parameters\n", - "obs = {'instrument': 'WFI',\n", - " 'filters': ['F129'],\n", - " 'detectors': 1,\n", - " 'background': 'pandeia',\n", - " 'observations_id': 42,\n", - " 'exptime': 300,\n", - " 'offsets': [{'offset_id': 1,\n", - " 'offset_centre': False,\n", - " 'offset_ra': 0.0,\n", - " 'offset_dec': 0.0,\n", - " 'offset_pa': 0.0\n", - " }]\n", - " }" + "Feed the dictionary to an instance of the `ObservationModule` class, while also specifying the central coordinates (90 degrees right ascension and 30 degrees declination) and position angle (0 degrees) as keyword arguments." ] }, { @@ -142,12 +133,31 @@ "metadata": {}, "outputs": [], "source": [ - "# Build observation parameters\n", - "obs = {'instrument': 'WFI', 'filters': ['F129'], 'detectors': 1,\n", - " 'background': 'pandeia', 'observations_id': 42, 'exptime': 300,\n", - " 'offsets': [{'offset_id': 1, 'offset_centre': False, 'offset_ra': 0.0,\n", - " 'offset_dec': 0.0, 'offset_pa': 0.0}]\n", - " }" + "obm = ObservationModule(obs, ra=90, dec=30, pa=0, seed=42, cores=6)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Create a simple astronomical scene\n", + "\n", + "The other requirement is the FITS source catalog. STIPS accepts [several types of tables and catalogs](https://stips.readthedocs.io/en/latest/using_stips/catalogue_formats.html#stips-table-formats) in either the IPAC text or FITS BinTable format. This example uses a Mixed Catalog, which requires the following columns:\n", + "\n", + "* `id`: Object ID\n", + "* `ra`: Right ascension (RA), in degrees\n", + "* `dec`: Declination (DEC), in degrees\n", + "* `flux`: Flux, in `units` (defined below)\n", + "* `type`: Approximation used to profile a source. Options are `'sersic'` (for extended sources) or `point`\n", + "* `n`: Sersic profile index\n", + "* `phi`: Position angle of the Sersic profle's major axis, in degrees\n", + "* `ratio`: Ratio of the Sersic profile's major and minor axes\n", + " * Since `n`, `phi`, `ratio` only apply to extended sources, they are ignored in rows where `type` is set to `'point'`.\n", + "* `notes`: Optional per-source comments\n", + "* `units`: One of `‘p’` for photons/s, `‘e’` for electrons/s, `‘j’` for Jansky, or `‘c’` for counts/s\n", + "\n", + "The catalog we create below contains two sources located near the central coordinates specified in the `ObservationModule`." ] }, { @@ -156,61 +166,63 @@ "metadata": {}, "outputs": [], "source": [ - "# Create observation object\n", - "obm = ObservationModule(obs, ra=90, dec=30, pa=0, seed=42, cores=6)" + "cols = [\n", + " fits.Column(name='id', array=[1, 2], format='K'),\n", + " fits.Column(name='ra', array=[90.02, 90.03], format='E'),\n", + " fits.Column(name='dec', array=[29.98, 29.97], format='D'),\n", + " fits.Column(name='flux', array=[0.00023, 0.0004], format='D'),\n", + " fits.Column(name='type', array=['point', 'point'], format='8A'),\n", + " fits.Column(name='n', array=[0, 0], format='D'),\n", + " fits.Column(name='re', array=[0, 0], format='D'),\n", + " fits.Column(name='phi', array=[0, 0], format='D'),\n", + " fits.Column(name='ratio', array=[0, 0], format='D'),\n", + " fits.Column(name='notes', array=['', ''], format='8A'),\n", + " fits.Column(name='units', array=['j', 'j'], format='8A')\n", + "]" ] }, { - "cell_type": "code", - "execution_count": 1, + "cell_type": "markdown", "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "'################################################################################'" - ] - }, - "execution_count": 1, - "metadata": {}, - "output_type": "execute_result" - } - ], "source": [ - "'#' * 80" + "Next, use the columns created above to create an output FITS table in the BinTable format and assign header keys that specify the filter and catalog type mentioned earlier." ] }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "rm ../webbpsf_tutorial/webbpsf_tutorial.ipynb ../pandeia_tutorial/pandeia_tutorial.ipynb" + "hdut = fits.BinTableHDU.from_columns(cols)\n", + "hdut.header['TYPE'] = 'mixed'\n", + "hdut.header['FILTER'] = 'F129'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "#### Create a simple astronomical scene" + "From there, save the table to disk." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cat_file = 'catalog.fits'\n", + "hdut.writeto(cat_file, overwrite=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "* `id`: Object ID\n", - "* `ra`: Right ascension (RA), in degrees\n", - "* `dec`: Declination (DEC), in degrees\n", - "* `flux`: Flux, in `units` (defined below)\n", - "* `type`: Type of PROFILE(???) (`'point'` or `'sersic'`)\n", - "* `n`: Sersic profile index (in what system???)\n", - "* `re`: Half-light radius, in pixels\n", - "* `phi`: Position angle of major axis, in degrees\n", - "* `ratio`: Major and minor axis ratio\n", - "* `notes`: Per-source comments\n", - "* `units`: ???" + "#### Simulate an image\n", + "\n", + "With the observation module and source catalog in tow, STIPS can take over the image simulation process. First, trigger the initialization of a new observation." ] }, { @@ -219,36 +231,43 @@ "metadata": {}, "outputs": [], "source": [ - "from astropy.io import fits\n", - " \n", - "cols = []\n", - "cols.append(fits.Column(name='id', array=[1, 2], format='K'))\n", - "cols.append(fits.Column(name='ra', array=[90.02, 90.03], format='E'))\n", - "cols.append(fits.Column(name='dec', array=[29.98, 29.97], format='D'))\n", - "cols.append(fits.Column(name='flux', array=[0.00023, 0.0004], format='D'))\n", - "cols.append(fits.Column(name='type', array=['point', 'point'], format='8A'))\n", - "cols.append(fits.Column(name='n', array=[0, 0], format='D'))\n", - "cols.append(fits.Column(name='re', array=[0, 0], format='D'))\n", - "cols.append(fits.Column(name='phi', array=[0, 0], format='D'))\n", - "cols.append(fits.Column(name='ratio', array=[0, 0], format='D'))\n", - "cols.append(fits.Column(name='notes', array=['', ''], format='8A'))\n", - "cols.append(fits.Column(name='units', array=['j', 'j'], format='8A'))\n", - " \n", - "# Create output fits table\n", - "hdut = fits.BinTableHDU.from_columns(cols)\n", - "hdut.header['TYPE'] = 'mixed'\n", - "hdut.header['FILTER'] = 'F129'\n", - " \n", - "# Write to disk\n", - "cat_file = 'catalog.fits'\n", - "hdut.writeto(cat_file, overwrite=True)" + "obm.nextObservation()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, simulate an image containing sources from the catalog you created earlier. Then, and add error residuals to the image." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "cat_name = obm.addCatalogue(cat_file)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "obm.addError()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "#### Simulate an image" + "Finish by saving the outputs to a FITS file. `ObservationModule`'s `finalize()` method can also return a mosaic and a list of the simulation's initial parameters." ] }, { @@ -257,14 +276,8 @@ "metadata": {}, "outputs": [], "source": [ - "# Initialize the local instrument\n", - "obm.nextObservation()\n", - " \n", - "# Add catalog with sources\n", - "cat_name = obm.addCatalogue(cat_file)\n", - " \n", - "# Add error to image\n", - "psf_file = obm.addError()" + "fits_file_1, _, params_1 = obm.finalize(mosaic=False)\n", + "print(f\"Output FITS file is {fits_file_1}\")" ] }, { @@ -273,9 +286,8 @@ "metadata": {}, "outputs": [], "source": [ - "# Call the final method\n", - "fits_file, mosaic_file, params = obm.finalize(mosaic=False)\n", - "print(f\"Output FITS file is {fits_file}\")" + "for prm in params_1:\n", + " print(prm)" ] }, { @@ -284,14 +296,27 @@ "metadata": {}, "outputs": [], "source": [ - "# DISPLAY IMAGE" + "img_1 = fits.getdata(fits_file_1)[1000:1500, 2500:3000]\n", + "\n", + "plt.imshow(img_1, vmax=1.5, origin='lower', cmap='bone')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Generate scenes from user-created catalogs\n", + "\n", + "STIPS can simulate scenes by importing pre-existing catalogs (as in the first example) or by using built-in functionality that generates collections of stars or galaxies based on user-specified parameters. Below, we will use the STIPS `SceneModule` to generate a stellar population and a galactic population." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Generate scene from a user-created catalog" + "#### Create a stellar population\n", + "\n", + "First, we specify parameters that will be the same in both scenes and use them to initalize separate `SceneModule` instances." ] }, { @@ -300,110 +325,218 @@ "metadata": {}, "outputs": [], "source": [ - "# USE SAME OBS PARAMS FOR BOTH CATALOGS\n", - "obs_prefix = 'notebook_example'\n", + "obs_prefix_1 = 'notebook_example1'\n", "obs_ra = 150.0\n", - "obs_dec = -2.5" + "obs_dec = -2.5\n", + "\n", + "scm_stellar = SceneModule(out_prefix=obs_prefix_1, ra=obs_ra, dec=obs_dec)\n", + "scm_galactic = SceneModule(out_prefix=obs_prefix_1, ra=obs_ra, dec=obs_dec)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, create a dictionary containing parameters of the desired stellar population to pass to one of the `SceneModule` instances' `CreatePopulation()` method." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "stellar_parameters = {'n_stars': 100, 'age_low': 7.5e12, 'age_high': 7.5e12,\n", + " 'z_low': -2.0, 'z_high': -2.0, 'imf': 'salpeter',\n", + " 'alpha': -2.35, 'binary_fraction': 0.1,\n", + " 'distribution': 'invpow', 'clustered': True,\n", + " 'radius': 100.0, 'radius_units': 'pc',\n", + " 'distance_low': 20.0, 'distance_high': 20.0,\n", + " 'offset_ra': 0.0, 'offset_dec': 0.0\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "_(Find a full accounting of each dictionary entry's meaning in the `CreatePopulation()` docstring.)_" ] }, { "cell_type": "code", "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "scm_stellar.CreatePopulation?" + ] + }, + { + "cell_type": "markdown", "metadata": {}, + "source": [ + "Pass the dictionary to the proper `SceneModule` instance's `CreatePopulation()` method. Running the method will save the newly-generated population to disk." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ - "# STELLAR\n", - "scm_stellar = SceneModule(out_prefix=obs_prefix, ra=obs_ra, dec=obs_dec)\n", - " \n", - "stellar_parameters = {\n", - " 'n_stars': 100,\n", - " 'age_low': 7.5e12,\n", - " 'age_high': 7.5e12,\n", - " 'z_low': -2.0,\n", - " 'z_high': -2.0,\n", - " 'imf': 'salpeter',\n", - " 'alpha': -2.35,\n", - " 'binary_fraction': 0.1,\n", - " 'clustered': True,\n", - " 'distribution': 'invpow',\n", - " 'radius': 100.0,\n", - " 'radius_units': 'pc',\n", - " 'distance_low': 20.0,\n", - " 'distance_high': 20.0,\n", - " 'offset_ra': 0.0,\n", - " 'offset_dec': 0.0\n", - " }\n", - " \n", "stellar_cat_file = scm_stellar.CreatePopulation(stellar_parameters)\n", "print(f\"Stellar population saved to file {stellar_cat_file}\")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Create a galactic population\n", + "\n", + "Repeat the population generation process, now by creating a dictionary containing parameters of a desired _galactic_ population, passing it to the other `SceneModule` instance's `CreateGalaxies()` method, and saving that result to disk." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "galaxy_parameters = {'n_gals': 10, 'z_low': 0.0, 'z_high': 0.2,\n", + " 'rad_low': 0.01, 'rad_high': 2.0,\n", + " 'sb_v_low': 30.0, 'sb_v_high': 25.0,\n", + " 'distribution': 'uniform', 'clustered': False,\n", + " 'radius': 200.0, 'radius_units': 'arcsec',\n", + " 'offset_ra': 0.0, 'offset_dec': 0.0\n", + " }" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "scm_galactic.CreateGalaxies??" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# GALACTIC\n", - "scm_galactic = SceneModule(out_prefix=obs_prefix, ra=obs_ra, dec=obs_dec)\n", - " \n", - "galaxy_parameters = {\n", - " 'n_gals': 10,\n", - " 'z_low': 0.0,\n", - " 'z_high': 0.2,\n", - " 'rad_low': 0.01,\n", - " 'rad_high': 2.0,\n", - " 'sb_v_low': 30.0,\n", - " 'sb_v_high': 25.0,\n", - " 'distribution': 'uniform',\n", - " 'clustered': False,\n", - " 'radius': 200.0,\n", - " 'radius_units': 'arcsec',\n", - " 'offset_ra': 0.0,\n", - " 'offset_dec': 0.0,\n", - " }\n", - " \n", "galaxy_cat_file = scm_galactic.CreateGalaxies(galaxy_parameters)\n", "print(f\"Galaxy population saved to file {galaxy_cat_file}\")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Set up an observation (first pointing)\n", + "\n", + "Once you've created a scene, you can use STIPS to simulate as many exposures of it in as many orientations as you'd like. In STIPS, a single telescope pointing is called an _offset_, and a collection of exposures is an _observation_." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We start this subsection by creating a single _offset_ that is dithered by 2 degrees in right ascension and rotated by 0.5 degrees in position angle from the center of the scene." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "offset_1 = {'offset_id': 1, 'offset_centre': False,\n", + " # True would center each detector on the same on-sky point\n", + " 'offset_ra': 2.0, 'offset_dec': 0.0, 'offset_pa': 0.5\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "That offset is contained within an _observation_ that is taken with the F129 filter, uses detectors WFI01 through WFI03, and has an exposure time of 1500 seconds. We also apply distortion and specify a sky background of 0.24 counts/s/pixel." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "observation_parameters_1 = {'instrument': 'WFI', 'filters': ['F129'],\n", + " 'detectors': 3, 'distortion': True,\n", + " 'background': 0.24, 'observations_id': 1,\n", + " 'exptime': 1500, 'offsets': [offset_1]\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "STIPS can also apply various types of error residuals to the observation. Here, we only include residuals from flatfielding and the dark current." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "residuals_1 = {'residual_flat': True, 'residual_dark': True,\n", + " 'residual_cosmic': False, 'residual_poisson': False,\n", + " 'residual_readnoise': False\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, feed the dictionaries of observation and residual parameters to an instance of the `ObservationModule` class, alongside the observation prefix, right ascension, and declination specified earlier in this example." + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# HOW DOES THIS RELATE TO OTHERS?\n", - "offset_1 = {\n", - " 'offset_id': 1,\n", - " 'offset_centre': False,\n", - " 'offset_ra': 2.0,\n", - " 'offset_dec': 0.0,\n", - " 'offset_pa': 0.5\n", - " }\n", - " \n", - "residuals_1 = {\n", - " 'residual_flat': True,\n", - " 'residual_dark': True,\n", - " 'residual_cosmic': False,\n", - " 'residual_poisson': False,\n", - " 'residual_readnoise': False\n", - " }\n", - " \n", - "observation_parameters_1 = {\n", - " 'instrument': 'WFI',\n", - " 'filters': ['F129'],\n", - " 'detectors': 3,\n", - " 'distortion': True,\n", - " 'background': 0.24,\n", - " 'observations_id': 1,\n", - " 'exptime': 1500,\n", - " 'offsets': [offset]\n", - " }\n", - " \n", - "obm_1 = ObservationModule(observation_parameters_1, out_prefix=obs_prefix,\n", - " ra=obs_ra, dec=obs_dec, residual=residuals_1)\n", - " \n", + "obm_1 = ObservationModule(observation_parameters_1, residual=residuals_1,\n", + " out_prefix=obs_prefix_1, ra=obs_ra, dec=obs_dec)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Call the `ObservationModule` object's `nextObservation()` method to move to the first offset/filter combination (`offset_1` and F129) contained in the object." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ "obm_1.nextObservation()" ] }, @@ -411,7 +544,135 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Modify observation parameters" + "Now that the observation is ready, add the stellar and galactic populations to it. (Each population may take about a minute to load.) Follow up by adding in the sources of error specified upon its initialization." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "output_stellar_catalogs_1 = obm_1.addCatalogue(stellar_cat_file)\n", + "output_galaxy_catalogs_1 = obm_1.addCatalogue(galaxy_cat_file)\n", + "print(f\"Output Catalogs are {output_stellar_catalogs_1} and \"\n", + " f\"{output_galaxy_catalogs_1}.\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "obm_1.addError()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As before, finish by saving the simulated image to a FITS file." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fits_file_1, _, params_1 = obm_1.finalize(mosaic=False)\n", + "print(f\"Output FITS file is {fits_file_1}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "img_1 = fits.getdata(fits_file_1, ext=1)\n", + "\n", + "norm_1 = mpl.colors.LogNorm(vmin=img_1.min(), vmax=img_1.max())\n", + "plt.imshow(img_1, norm=norm_1, origin='lower', cmap='bone')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Modify an observation (second pointing)\n", + "\n", + "To observe the same scene under different conditions, make a new `ObservationModule` object that takes updated versions of the input observation parameter and residual dictionaries. Taken collectively, the collection of resulting `ObservationModule` objects can be thought of as a dithered set of observations." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For the second observation, create another offset that is dithered by 10 degrees in right ascension and rotated by 27 degrees in position angle from the center of the scene." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "offset_2 = {'offset_id': 1, 'offset_centre': False,\n", + " # True centers each detector on same point\n", + " 'offset_ra': 10.0, 'offset_dec': 0.0, 'offset_pa': 27\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The second observation is the same as the first (F129 filter, detectors WFI01 through WFI03, exposure time of 1500 seconds, distortion, and a sky background of 0.24 counts/s/pixel) besides the substitution of the first set of offset parameters for the second." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "observation_parameters_2 = {'instrument': 'WFI', 'filters': ['F129'],\n", + " 'detectors': 3, 'distortion': True,\n", + " 'background': 0.24, 'observations_id': 1,\n", + " 'exptime': 1500, 'offsets': [offset_2]\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This time, we include residuals from the flatfield (like before) and read noise (unlike before)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "residuals_2 = {'residual_flat': True, 'residual_dark': False,\n", + " 'residual_cosmic': False, 'residual_poisson': False,\n", + " 'residual_readnoise': True\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create the new `ObservationModule` object and initialize it for a simulation." ] }, { @@ -420,35 +681,19 @@ "metadata": {}, "outputs": [], "source": [ - "offset_2 = {\n", - " 'offset_id': 1,\n", - " 'offset_centre': False,\n", - " 'offset_ra': 10.0,\n", - " 'offset_dec': 0.0,\n", - " 'offset_pa': 27\n", - " }\n", - " \n", - "residuals_2 = {\n", - " 'residual_flat': True,\n", - " 'residual_dark': False,\n", - " 'residual_cosmic': False,\n", - " 'residual_poisson': False,\n", - " 'residual_readnoise': True\n", - " }\n", - " \n", - "observation_parameters_2 = {\n", - " 'instrument': 'WFI',\n", - " 'filters': ['F129'],\n", - " 'detectors': 3,\n", - " 'distortion': True,\n", - " 'background': 0.24,\n", - " 'observations_id': 1,\n", - " 'exptime': 1500,\n", - " 'offsets': [offset]\n", - " }\n", - " \n", - "obm_2 = ObservationModule(observation_parameters_2, out_prefix=obs_prefix, ra=obs_ra, dec=obs_dec, residual=residuals_2)\n", - " \n", + "obs_prefix_2 = 'notebook_example2'\n", + "obm_2 = ObservationModule(observation_parameters_2, residuals=residuals_2,\n", + " out_prefix=obs_prefix_2, ra=obs_ra, dec=obs_dec)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ "obm_2.nextObservation()" ] }, @@ -456,7 +701,81 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Add artifical point source to an observation" + "Add the stellar and galactic populations to the new observation along with the sources of error chosen above." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "output_stellar_catalogs_2 = obm_2.addCatalogue(stellar_cat_file)\n", + "output_galaxy_catalogs_2 = obm_2.addCatalogue(galaxy_cat_file)\n", + "print(f\"Output Catalogs are {output_stellar_catalogs_2} and \"\n", + " f\"{output_galaxy_catalogs_2}.\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "obm_2.addError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fits_file_2, _, params_2 = obm_2.finalize(mosaic=False)\n", + "print(f\"Output FITS file is {fits_file_2}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Below, see the two pointings side by side." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "img_2 = fits.getdata(fits_file_2, ext=1)\n", + "norm_2 = mpl.colors.LogNorm(vmin=img_2.min(), vmax=img_2.max())\n", + "\n", + "fig_2, ax_2 = plt.subplots(1, 2, figsize=(8,4))\n", + "ax_2[0].imshow(img_1, norm=norm_1, origin='lower', cmap='bone')\n", + "ax_2[1].imshow(img_2, norm=norm_2, origin='lower', cmap='bone')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Add artifical point source to an observation\n", + "\n", + "With the STIPS `makePSF` utility, users can \"clip\" a PSF from a given detector pixel position in a scene and inject it elsewhere.\n", + "\n", + "This is achieved by 1) applying bi-linear interpolation of the 3x3 PSF library array to compute the best PSF at the specified integer SCA pixels, then 2) performing bicubic interpolations over the PSF's supersampled pixel grid to fill out its sub-pixel positions. The resulting PSF can then be injected elsewhere in an existing scene or used to create new scenes." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The library PSF is made from detector SCA01 with the F129 filter. It is in the same directory as this notebook, so we open it below." ] }, { @@ -466,7 +785,40 @@ "outputs": [], "source": [ "with fits.open('psf_WFI_2.0.0_F129_sca01.fits') as hdul:\n", - " test_psf = stips.utilities.makePSF.make_epsf(hdul[0].data)" + " test_psf = stips.utilities.makePSF.make_epsf(hdul[0].data[0])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.imshow(test_psf, cmap='bone')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Specify the pixel coordinates of the source's center (which are just the center of the PSF in this case) and the pixel length of the intended cutout (1/8th of the full example PSF's pixel length)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "psf_mid_pixel = (test_psf.shape[0] - 1) // 2\n", + "boxsize = test_psf.shape[0] // 8" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Inject the source at pixel (2000, 2000) with a flux of 3000 count/s _**(??)**_ in a copy of the simulated image from the previous section's second (final) observation. Can you see the newly injected PSF in the comparison plot?" ] }, { @@ -475,8 +827,14 @@ "metadata": {}, "outputs": [], "source": [ - "# VARS SEEM DIFFERENT – DOES THIS NEED CODE FROM PREVIOUS EXAMPLES? IF NOT, RDOX REFERENCE\n", - "# TO COPYING EXAMPLE 2 CODE IS WRONG." + "img_3 = img_2.copy()\n", + "xpix = 2000\n", + "ypix = 2000\n", + "flux = 3000\n", + "\n", + "img_3_inj = stips.utilities.makePSF.place_source(xpix, ypix, flux, img_3,\n", + " test_psf, boxsize=boxsize,\n", + " psf_center=psf_mid_pixel)" ] }, { @@ -485,45 +843,14 @@ "metadata": {}, "outputs": [], "source": [ - "offset = {\n", - " 'offset_id': 1,\n", - " 'offset_centre': False,\n", - " 'offset_ra': 0.0,\n", - " 'offset_dec': 0.0,\n", - " 'offset_pa': 0.0\n", - " }\n", - " \n", - "observation_parameters = {\n", - " 'instrument': 'WFI',\n", - " 'filters': ['F129'],\n", - " 'detectors': 1,\n", - " 'distortion': False,\n", - " 'background': 0.15,\n", - " 'observations_id': 1,\n", - " 'exptime': 1000,\n", - " 'offsets': [offset]\n", - " }\n", - " \n", - "obm = ObservationModule(observation_parameters, out_prefix=obs_prefix, ra=obs_ra, dec=obs_dec)\n", - " \n", - "# nextObservation is called to move between different combinations of offset and filter.\n", - "# It must be called once in order to initialize the observation module to the first observation before adding catalogs.\n", - "obm.nextObservation()\n", - " \n", - "output_stellar_catalogs = obm.addCatalogue(stellar_cat_file)\n", - "output_galaxy_catalogs = obm.addCatalogue(galaxy_cat_file)\n", - " \n", - "print(f\"Output Catalogs are {output_stellar_catalogs} and {output_galaxy_catalogs}.\")\n", - " \n", - "psf_file = obm.addError()\n", - " \n", - "print(\"PSF File is {psf_file}\")\n", - " \n", - "fits_file, mosaic_file, params = obm.finalize(mosaic=False)\n", - " \n", - "print(f\"Output FITS file is {fits_file}\")\n", - "print(f\"Output Mosaic File is {mosaic_file}\")\n", - "print(f\"Observation Parameters are {params}\")" + "norm_3 = mpl.colors.LogNorm(vmin=img_3_inj.min(), vmax=img_3_inj.max())\n", + "\n", + "fig_3, ax_3 = plt.subplots(1, 2, figsize=(8,4))\n", + "ax_3[0].imshow(img_2, norm=norm_2, origin='lower', cmap='bone')\n", + "ax_3[1].imshow(img_3_inj, norm=norm_3, origin='lower', cmap='bone')\n", + "\n", + "for ax in ax_3:\n", + " ax.add_patch(plt.Circle((xpix, ypix), 75, color='r', alpha=.7, fill=False))" ] }, { @@ -586,7 +913,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.5" + "version": "3.11.9" } }, "nbformat": 4,