diff --git a/notebooks/concepts/generate_mosviz_photometry.ipynb b/notebooks/concepts/generate_mosviz_photometry.ipynb index 53037a6321..a30c6386d0 100644 --- a/notebooks/concepts/generate_mosviz_photometry.ipynb +++ b/notebooks/concepts/generate_mosviz_photometry.ipynb @@ -2,24 +2,39 @@ "cells": [ { "cell_type": "markdown", - "id": "e5e8735c-0e9a-4533-a4e2-f2b6b9ece9f9", + "id": "01609adc-98ac-41d2-8030-8eb545ca2fb0", "metadata": {}, "source": [ "# Synthetic image creation for MOSviz pipeline data\n", "\n", - "**Motiviation**: The synthetic dataset we've received from the JWST data pipeline team for use in MOSviz contains simulated spectra from NIRSpec, but no simulated photometry from NIRCam. We'd like to have test imagery to display in MOSviz alongside the 2D and 1D spectra we already possess. We've learned that the pipeline team has no plans to produce any, so we attempt to piece together our own.\n", + "**Motiviation**: The synthetic dataset we currently possess from the JWST data pipeline team for use in MOSviz contains simulated spectra from NIRSpec but no associated NIRCam photometry. We would like to have test imagery to display in MOSviz alongside these 2D and 1D spectra. We are unsure whether the pipeline team has plans to produce any.\n", "\n", - "**Goal**: Populate an image of background noise with properly-scaled galaxy cutouts sourced from an HST image and placed at their respective WCS locations. These galaxies' real spectra don't necessarily correspond with those in our dataset, but at this point we care more about the veneer of having photometry to match with our spectra.\n", + "**Goal**: Populate an image of background noise with properly-scaled galaxy cutouts sourced from a Hubble Space Telescope image and placed at their analogous locations in the new image. These galaxies' real spectra do not necessarily correspond with those in our dataset, but we care more about the veneer of having photometry to match with our spectra at this point.\n", "\n", "**Execution**: We pull our galaxy cutouts and catalog information from the Hubble Deep Field image. ~~ASTRODEEP's image of the [Abell 2774 Parallel](http://astrodeep.u-strasbg.fr/ff/?img=JH140?cm=grayscale) | [MACS J0416.1-2403 Parallel](http://astrodeep.u-strasbg.fr/ff/?ffid=FF_M0416PAR&id=1264&cm=grayscale)~~. We sought to use [Artifactory](https://bytesalad.stsci.edu/ui/repos/tree/General/jwst-pipeline%2Fdev%2Ftruth) to obtain a range of RA/Dec over which to project our synthetic image, but that information was absent. Instead, we place the image over a manually chosen RA/Dec range and place the cutouts in randomly selected locations within that field of view.\n", "\n", "**Issues**:\n", - "- We wanted to scrape RA/Dec information the data pipeline products to get a range of coordinates over which to scale our synthetic image, but it appears that the pipeline's data products don't have `\"TARG_RA\"` or `\"TARG_DEC\"` keywords in their headers.\n", + "- We wanted to scrape RA/Dec information the data pipeline products to get a range of coordinates over which to scale our synthetic image, but it appears that the pipeline's data products lack `\"TARG_RA\"` or `\"TARG_DEC\"` keywords in their headers.\n", " - The data products also don't appear to have WCS information. We don't strictly need it to achieve this notebook's goals, but it would be convenient to have.\n", - "- There doesn't appear to be an observation with level 2 data in `jwst-pipeline/truth/test_nirspec_mos_spec2` and level 3 data in `jwst-pipeline/truth/test_nirspec_mos_spec3`. All observations are either level 2 only or level 3 only.\n", - "- _(Resolved)_ A good number of the cutouts from the first couple of field images we tested had intrustions from other galaxies due to crowding. We settled on the Hubble Deep Field as a good source image, but if we hadn't, we may have considered using galaxies modeled with Sersic profiles to get cleaner cutouts to inject into our synthetic image.\n", + "- There does not appear to be an observation with level 2 data in `jwst-pipeline/truth/test_nirspec_mos_spec2` and level 3 data in `jwst-pipeline/truth/test_nirspec_mos_spec3`. All observations are either level 2 only or level 3 only.\n", + "- _(Resolved)_ A good number of the cutouts from the first couple of field images we tested had intrustions from other galaxies due to crowding. We settled on the Hubble Deep Field as a good source image, but had we not, we may have considered using galaxies modeled with Sersic profiles to get cleaner cutouts to inject into our synthetic image." + ] + }, + { + "cell_type": "markdown", + "id": "57df882f-e15d-4f16-931c-9f495893c0e1", + "metadata": {}, + "source": [ + "### Import packages\n", "\n", - "*Writers: Robel Geda and O. Justin Otor*" + "- We use `astropy.io.fits` to read in existing FITS files and write a new one with the synthetic image.\n", + "- The objects from `astropy.nddata` help with creating cutouts once we've identified galaxies we'd like to take from the field image.\n", + "- The methods from `astropy.stats` work with image data that's clipped to within a certain number of deviations from the mean.\n", + "- The objects from `astropy.table` help with reading an modifiying tabular data.\n", + "- `astropy.wcs.WCS` creates a World Coordinate System object that's useful for transforming on-sky data from one patch to another.\n", + "- `glob.glob` lists local files that match a given pattern.\n", + "- We use `matplotlib.pyplot` to preview the field image, the cutouts, and various stages of our synthetic image.\n", + "- We use `numpy` to facilitate several specialized mathematical and array-related operations." ] }, { @@ -45,7 +60,9 @@ "id": "fffd44de-ee89-4a46-8c75-198a94cf124d", "metadata": {}, "source": [ - "### Generate galaxy cutouts" + "### Generate galaxy cutouts\n", + "\n", + "The galaxy cutouts come from the Hubble Deep Field image. To retrieve them, we download the image itself and its associated catalog data, search the catalog for the brightest galaxies, then find those galaxies in the image." ] }, { @@ -70,15 +87,17 @@ "metadata": {}, "outputs": [], "source": [ - "# download those sources' locations in the image and sort them by brightness\n", + "# download sources' location and flux information\n", "source_info1 = Table.read('https://archive.stsci.edu/pub/hlsp/hdf/wfpc_hdfn_v2catalog/HDFN_wfpc_v2generic.cat',\n", " format='ascii')\n", "source_info2 = Table.read('https://archive.stsci.edu/pub/hlsp/hdf/wfpc_hdfn_v2catalog/HDFN_f814_v2.cat',\n", " format='ascii')\n", - "sources = table_join(source_info1, source_info2)\n", + "sources = join(source_info1, source_info2)\n", "\n", - "# confirm that both tables contain the same objects in the same order\n", - "(source_info1['NUMBER'] == source_info2['NUMBER']).sum() == len(source_info1) == len(source_info2)" + "# confirm that both tables contain the same objects in the same order (True)\n", + "( (source_info1['NUMBER'] == source_info2['NUMBER']).sum()\n", + " == len(source_info1)\n", + " == len(source_info2) )" ] }, { @@ -88,10 +107,8 @@ "metadata": {}, "outputs": [], "source": [ - "# sort sources by flux contained within 71.1 pixel diameter of source (11)\n", + "# sort sources by flux within 71.1 pixel diameter of source, or aperture 11\n", "sources.sort('FLUX_APER_11', reverse=True)\n", - "#sources.sort('KRON_RADIUS', reverse=True)\n", - "#sources.sort('FWHM_IMAGE')\n", "\n", "# filter out likely stars and sources with negative flux\n", "sources = sources[(sources['CLASS_STAR'] < .5)\n", @@ -121,6 +138,14 @@ " sources['DELTA_J2000'])" ] }, + { + "cell_type": "markdown", + "id": "6b4f1421-e3c3-4112-a475-dc1b7ad07146", + "metadata": {}, + "source": [ + "Note: Depending on the value of `catalog_size`, the following cell can produce a lot of output. Right-click the cell and select \"Enable Scrolling for Outputs\" to expand it or \"Disable Scrolling for Outputs\" to condense it." + ] + }, { "cell_type": "code", "execution_count": null, @@ -162,6 +187,14 @@ " break" ] }, + { + "cell_type": "markdown", + "id": "c3504d50-99cf-4d73-9678-1eff2ca43676", + "metadata": {}, + "source": [ + "We also save image statistics calculated from pixels within a chosen number of standard deviations from the image's mean intensity. Some of them may be useful in creating the synthetic image later on." + ] + }, { "cell_type": "code", "execution_count": null, @@ -169,7 +202,8 @@ "metadata": {}, "outputs": [], "source": [ - "clipped_mean, clipped_median, clipped_stddev = sigma_clipped_stats(image_data, sigma=3.)" + "clipped_mean, clipped_median, clipped_stddev = sigma_clipped_stats(image_data,\n", + " sigma=3.)" ] }, { @@ -177,7 +211,11 @@ "id": "53891064-621c-4f6f-91da-3a24741ad632", "metadata": {}, "source": [ - "### Extract destination RA/Dec from spectra files" + "### Extract destination RA/Dec from spectra files\n", + "\n", + "Note that the following cells only run if on a connection with access to the STScI VPN. They use copies of JWST pipeline files from May 19, 2020.\n", + "\n", + "To run them elsewhere and/or use more current files, visit the [Artifactory](https://bytesalad.stsci.edu/ui/repos/tree/General/jwst-pipeline%2Fdev%2Ftruth%2Ftest_nirspec_mos_spec3), ensure the `jwst-pipeline/dev/truth/test_nirspec_mos_spec3/` folder is selected in the directory tree, and download all files whose names begin with `jw00626-o030` and end with `nirspec_f170lp-g235m_x1d.fits`. (The original version of the notebook uses six of these Level 3 data products.) Then, change `filepath` in the next cell to the location where you saved the downloaded files. Remember the trailing forward slash." ] }, { @@ -187,45 +225,38 @@ "metadata": {}, "outputs": [], "source": [ - "# save level 3 spectra FITS header information\n", - "# (is there a way to do so using the URL?)\n", - "# 'https://bytesalad.stsci.edu/ui/repos/tree/General/jwst-pipeline%2Fdev%2Ftruth%2Ftest_nirspec_mos_spec3%2Fjw00626-o030_s00000_nirspec_f170lp-g235m_s2d.fits'\n", - "x1d_header = fits.getheader('/Users/jotor/Downloads/jw00626-o030_s00000_nirspec_f170lp-g235m_x1d.fits')\n", - "s2d_header = fits.getheader('/Users/jotor/Downloads/jw00626-o030_s00000_nirspec_f170lp-g235m_s2d.fits')" + "# view level 3 spectra FITS header information\n", + "filepath = '/user/jotor/jwst-pipeline-lvl3/'\n", + "x1d_header = fits.getheader(filepath + 'jw00626-o030_s00000_nirspec_f170lp-g235m_x1d.fits')\n", + "#s2d_header = fits.getheader(filepath + 'jw00626-o030_s00000_nirspec_f170lp-g235m_s2d.fits')" ] }, { "cell_type": "code", "execution_count": null, - "id": "bb2858ac-e8e9-4c80-a2ac-72947ed7195b", + "id": "96dda8bd-5010-45ef-a937-b99ee54b8091", "metadata": {}, "outputs": [], "source": [ - "for k, v in x1d_header.items():\n", - " if 'or' in k.lower():\n", - " print(k, v)" + "x1d_header['TARG_RA'], x1d_header['TARG_DEC']" ] }, { "cell_type": "code", "execution_count": null, - "id": "96dda8bd-5010-45ef-a937-b99ee54b8091", + "id": "ade024d2-9788-4572-9304-af69c30e86da", "metadata": {}, "outputs": [], "source": [ - "(x1d_header['TARG_RA'], x1d_header['TARG_DEC'],\n", - " s2d_header['TARG_RA'], s2d_header['TARG_DEC'])" + "WCS(x1d_header)" ] }, { - "cell_type": "code", - "execution_count": null, - "id": "a4ecfbe4-bc85-4d94-b660-e50d0bf621ba", + "cell_type": "markdown", + "id": "809be77d-b099-40a3-8a85-42a7a2cf3244", "metadata": {}, - "outputs": [], "source": [ - "# no WCS data... not what we want\n", - "WCS(x1d_header)" + "Notice that this header lacks WCS information. Additionally, examining all of this observation's file headers reveals that they all have the same RA/Dec, which is not what we expect for its different \"pointings.\"" ] }, { @@ -236,13 +267,21 @@ "outputs": [], "source": [ "# search for RA/Dec information from Artifactory observation files\n", - "x1d_header_list = [fits.getheader(file) for file in glob('/Users/jotor/Downloads/jw00626*x1d.fits')]\n", + "x1d_header_list = [fits.getheader(file)\n", + " for file in glob(filepath + 'jw00626*x1d.fits')]\n", "\n", - "# all the same... not what we want\n", "ras, decs = np.array([[h['TARG_RA'], h['TARG_DEC']] for h in x1d_header_list]).T\n", "ras, decs" ] }, + { + "cell_type": "markdown", + "id": "debb3ccc-88ec-4e2d-bf50-69e8361f7d4b", + "metadata": {}, + "source": [ + "Since the headers lack the information we seek, we randomly generate our sources' RA/Dec information in predetermined patch of sky. We take the patch size from the size of NIRSpec Micro-Shutter Assembly (MSA) -- about 3.6'x3.4'-- approximated to a square field of view." + ] + }, { "cell_type": "code", "execution_count": null, @@ -250,8 +289,6 @@ "metadata": {}, "outputs": [], "source": [ - "# since all information is the same, we randomly generate our own RAs/decs\n", - "# (ranges based on NIRSpec MSA's on-sky projection size of 3.6x3.4 arcmins)\n", "np.random.seed(19)\n", "ras = np.random.uniform(0, 1/15, catalog_size)\n", "decs = np.random.uniform(-1/30, 1/30, catalog_size)" @@ -262,7 +299,9 @@ "id": "c4b5da62-c424-4603-adc3-4c550291d13a", "metadata": {}, "source": [ - "### Create synthetic image" + "### Create synthetic image\n", + "\n", + "We initialize a `numpy` array and fill it with normally-distributed background noise based on some of the clipped image statistics that were calculated earlier." ] }, { @@ -272,7 +311,6 @@ "metadata": {}, "outputs": [], "source": [ - "# create synthetic image onto which cutouts will be pasted\n", "synth_img_size = 1000\n", "synth_image = np.zeros((synth_img_size, synth_img_size))" ] @@ -286,9 +324,7 @@ "source": [ "# add noise\n", "synth_image += np.random.normal(loc=clipped_mean, scale=clipped_stddev*8,\n", - " size=synth_image.shape)\n", - "# synth_image += np.random.normal(loc=image_data.mean(), scale=image_data.std(),\n", - "# size=synth_image.shape)" + " size=synth_image.shape)" ] }, { @@ -298,7 +334,8 @@ "metadata": {}, "outputs": [], "source": [ - "plt.imshow(synth_image, cmap='bone')\n", + "imshow_params = {'cmap': 'bone', 'origin': 'lower'}\n", + "plt.imshow(synth_image, **imshow_params)\n", "plt.show()" ] }, @@ -307,7 +344,9 @@ "id": "674017bb-2cc1-437a-bf3a-449f306fa506", "metadata": {}, "source": [ - "### Fill out new WCS object for `synth_image`" + "### Fill out new WCS object for `synth_image`\n", + "\n", + "Creating a `WCS` object for `synth_image` allows it to be mathematically transformed into a projection on the sky. That projection can then be compared to other FITS images with their own `WCS` information." ] }, { @@ -322,27 +361,26 @@ ] }, { - "cell_type": "code", - "execution_count": null, - "id": "0d8aa070-1b6f-4a0a-9975-58bf5634ff86", - "metadata": { - "tags": [] - }, - "outputs": [], + "cell_type": "markdown", + "id": "c1a2876e-963b-4e62-8dd1-d68fb4b4e0a7", + "metadata": {}, "source": [ - "ra_bounds = np.array([ras.max(), ras.min()])\n", - "dec_bounds = np.array([decs.max(), decs.min()])" + "The next step is to calculate field of view information for `synth_image`." ] }, { "cell_type": "code", "execution_count": null, - "id": "80f9d82c-2a63-42f9-8cdc-2dc5d979062e", + "id": "068af407-9a8a-429c-9f7e-03f4f0a0b867", "metadata": {}, "outputs": [], "source": [ - "delta_ra = ras.max() - ras.min()\n", - "delta_dec = decs.max() - decs.min()" + "# find the range of sources in RA and dec\n", + "ra_bounds = np.array([ras.max(), ras.min()])\n", + "dec_bounds = np.array([decs.max(), decs.min()])\n", + "\n", + "delta_ra = np.ptp(ras)\n", + "delta_dec = np.ptp(decs)" ] }, { @@ -352,9 +390,9 @@ "metadata": {}, "outputs": [], "source": [ - "# set minimum FOV as maximum range in RA or dec\n", + "# save the maximum span in coordinates, RA or dec\n", "if delta_ra > delta_dec:\n", - " min_image_fov = abs(delta_ra * np.cos(np.pi/180 * dec_bounds.sum()/2))\n", + " min_image_fov = abs(delta_ra * np.cos(np.pi / 180 * dec_bounds.sum() / 2))\n", "else:\n", " min_image_fov = delta_dec\n", " \n", @@ -368,14 +406,22 @@ "metadata": {}, "outputs": [], "source": [ - "# scale this FOV by pixels\n", + "# scale this field of view (FOV) by pixels\n", "pix_scale = min_image_fov / synth_img_size\n", "\n", - "# add a buffer to the borders\n", + "# add a buffer to the FOV's borders to prevent clipping sources\n", "pix_scale *= 1.5\n", "pix_scale" ] }, + { + "cell_type": "markdown", + "id": "d28eb2e7-d3a1-4190-978a-8be1b30d3a3c", + "metadata": {}, + "source": [ + "With those calculations done, the `WCS` object is ready to be filled out." + ] + }, { "cell_type": "code", "execution_count": null, @@ -397,10 +443,20 @@ "synth_wcs" ] }, + { + "cell_type": "markdown", + "id": "5792df4f-c0d6-49a5-be72-a7c5ab39d9b4", + "metadata": {}, + "source": [ + "### Populate `synth_image` with the cutouts\n", + "\n", + "Finally, we convert the cutouts' coordinates to pixels and add them into `synth_image` to complete the creation of the mock image." + ] + }, { "cell_type": "code", "execution_count": null, - "id": "8f92545b-5c5d-4ab3-a9d6-ec61e2d329a6", + "id": "34c92d83-8e78-4cc5-8df3-728cb51d6567", "metadata": {}, "outputs": [], "source": [ @@ -409,14 +465,6 @@ "ras_pix, decs_pix" ] }, - { - "cell_type": "markdown", - "id": "5792df4f-c0d6-49a5-be72-a7c5ab39d9b4", - "metadata": {}, - "source": [ - "### Populate `synth_image` with the cutouts" - ] - }, { "cell_type": "code", "execution_count": null, @@ -439,11 +487,19 @@ "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(10,10))\n", - "ax.imshow(synth_image, vmin=0, vmax=synth_image.std()*3, origin='lower', cmap='bone')\n", - "#plt.imshow(synth_image, vmin=0, vmax=synth_image.mean()*3, origin='lower', cmap='bone')\n", + "ax.imshow(synth_image, vmin=0, vmax=synth_image.std()*3, **imshow_params)\n", + "\n", "plt.show()" ] }, + { + "cell_type": "markdown", + "id": "eb2a0045-74b6-4ac2-85a8-c535e1a64800", + "metadata": {}, + "source": [ + "Save the image, overwriting the previous file if it already exists in the current directory." + ] + }, { "cell_type": "code", "execution_count": null, @@ -454,6 +510,19 @@ "fits.writeto('synthetic_HDF_more.fits', synth_image,\n", " header=synth_wcs.to_header(), overwrite=True)" ] + }, + { + "cell_type": "markdown", + "id": "fe8817e2-51f4-4417-84ce-659c9449e87a", + "metadata": {}, + "source": [ + "

\n", + " Authors: Robel Geda and O. Justin Otor \n", + " \"Space\n", + "

\n", + "\n", + "-----" + ] } ], "metadata": {