From 98f0acabb68f3b83ae01dabb07e9df46b703ae18 Mon Sep 17 00:00:00 2001 From: Gagandeep Anand <30330521+gsanand@users.noreply.github.com> Date: Wed, 31 Jul 2024 17:14:50 -0400 Subject: [PATCH] Add files via upload --- .../align_multiple_visits.ipynb | 871 ++++++++++++++---- 1 file changed, 696 insertions(+), 175 deletions(-) diff --git a/notebooks/DrizzlePac/align_multiple_visits/align_multiple_visits.ipynb b/notebooks/DrizzlePac/align_multiple_visits/align_multiple_visits.ipynb index 22e1d7734..cecafea41 100644 --- a/notebooks/DrizzlePac/align_multiple_visits/align_multiple_visits.ipynb +++ b/notebooks/DrizzlePac/align_multiple_visits/align_multiple_visits.ipynb @@ -4,31 +4,75 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Optimizing Image Alignment for Multiple HST Visits" + "# Optimizing Image Alignment for Multiple HST Visits " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "
Note: The notebook in this repository 'Initialization.ipynb' goes over many of the basic concepts such as the setup of the environment/package installation and should be read first if you are new to HST images, DrizzlePac, or Astroquery.
" + "
This notebook requires creating and activating a virtual environment using the requirements file in this notebook's repository. Please also review the README file before using the notebook.
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Introduction\n", + "\n", + "## Table of Contents\n", "\n", - "This notebook demonstrates how to use `TweakReg` and `AstroDrizzle` tasks to align and combine images. While this example focuses on ACS/WFC data, the procedure can be almost identically applied to WFC3/UVIS images. This notebook is based on [a prior example](http://documents.stsci.edu/hst/HST_overview/documents/DrizzlePac/ch75.html) from the 2012 [DrizzlePac Handbook](http://documents.stsci.edu/hst/HST_overview/documents/DrizzlePac/toc.html), but has been updated for compatibility with the STScI AstroConda software distribution. There is a good deal of explanatory text in this notebook, and users new to DrizzlePac are encouraged to start with this tutorial. Additional DrizzlePac software documentation is available at [the readthedocs webpage](https://drizzlepac.readthedocs.io).\n", + "[Introduction](#intro)
\n", "\n", - "Before running the notebook, you will need to install the `astroquery` package, which is used to retrieve the data from the MAST archive and the `ccdproc` package which is used to query the image headers. \n", + "[1. Download the Observations with `astroquery`](#download)
\n", "\n", - "Summary of steps:\n", + "    [1.1 Check image header data](#check_keywords)
\n", + "    [1.2 Inspect the alignment](#check_wcs)
\n", "\n", - "1. Download the data from MAST using astroquery. \n", - "2. Align the images to a common reference frame, and update the WCS using `TweakReg`.\n", - "3. Combine the aligned images using `AstroDrizzle`." + "[2. Align with `TweakReg`](#tweakreg)
\n", + "\n", + "    [2.1 First test: Use 'default' parameters](#default)
\n", + "    [2.2 Second test: Set parameters for source finding](#sources)
\n", + "    [2.3 Inspect the quality of the fit](#fit_quality)
\n", + "\n", + "[3. Combine the Images using `AstroDrizzle`](#adriz)
\n", + "\n", + "    [3.1 Inspect the drizzled science image](#science)
\n", + "    [3.2 Inspect the final weight image](#wht)
\n", + "    [3.3 Discussion of `AstroDrizzle`](#discussion)
\n", + "\n", + "The following Python packages are required to run the Jupyter Notebook:\n", + " - [**glob**](https://docs.python.org/3/library/glob.html) - gather lists of filenames\n", + " - [**os**](https://docs.python.org/3/library/os.html) - change and make directories\n", + " - [**shutil**](https://docs.python.org/3/library/shutil.html) - perform high-level file operations\n", + "\n", + " - [**iPython**](https://ipython.readthedocs.io/en/stable/) - interactive handling\n", + " - [**matplotlib**](https://matplotlib.org) - create graphics\n", + " - [**numpy**](https://numpy.org) - math and array functions\n", + " \n", + " - [**astropy**](https://www.astropy.org) - FITS file handling and related operations\n", + " - [**astroquery**](https://astroquery.readthedocs.io/en/latest/) - download data and query databases\n", + " - [**drizzlepac**](https://drizzlepac.readthedocs.io/en/latest/) - align and combine HST images" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Introduction \n", + "[Table of Contents](#toc)\n", + "\n", + "This notebook demonstrates how to use `TweakReg` and `AstroDrizzle` tasks to align and combine images. While this example focuses on ACS/WFC data, the procedure can be almost identically applied to WFC3/UVIS images. This notebook is based on [a prior example](http://documents.stsci.edu/hst/HST_overview/documents/DrizzlePac/ch75.html) from the 2012 [DrizzlePac Handbook](http://documents.stsci.edu/hst/HST_overview/documents/DrizzlePac/toc.html), but has been updated for compatibility with the latest STScI software distribution. \n", + "\n", + "There is a good deal of explanatory text in this notebook, and users new to DrizzlePac are encouraged to start with this tutorial. Additional `DrizzlePac` software documentation is available at the [readthedocs webpage](https://drizzlepac.readthedocs.io).\n", + "\n", + "Before running the notebook, you will need to install the `astroquery` package, which is used to retrieve the data from the MAST archive. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Please also make sure to appropriately set your environment variables to allow Python to locate the appropriate reference files." ] }, { @@ -38,39 +82,77 @@ "outputs": [], "source": [ "import glob\n", - "import os\n", + "import os \n", "import shutil\n", "\n", - "from astropy.io import ascii\n", - "from astropy.io import fits\n", + "from collections import defaultdict\n", + "from IPython.display import clear_output\n", + "\n", + "import matplotlib.image as mpimg\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "from astropy.io import ascii, fits\n", "from astropy.table import Table\n", "from astropy.visualization import ZScaleInterval\n", + "\n", "from astroquery.mast import Observations\n", - "from drizzlepac import tweakreg\n", - "from drizzlepac import astrodrizzle\n", - "from IPython.display import Image\n", - "import matplotlib.pyplot as plt\n", "\n", - "%matplotlib inline" + "from drizzlepac import tweakreg, astrodrizzle\n", + "from drizzlepac.processInput import getMdriztabPars\n", + "\n", + "%matplotlib inline\n", + "%config InlineBackend.figure_format='retina'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## 1. Download the Data \n", + "Along with the above imports, we also set up the required local reference files." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Set the locations of reference files to your local computer.\n", + "os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'\n", + "os.environ['CRDS_SERVER'] = 'https://hst-crds.stsci.edu'\n", + "os.environ['CRDS_PATH'] = './crds_cache'\n", + "os.environ['iref'] = './crds_cache/references/hst/wfc3/'\n", + "os.environ['jref'] = './crds_cache/references/hst/acs/'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Download the Observations with `astroquery` \n", + "[Table of Contents](#toc)\n", + "\n", + "---\n", + "MAST queries may be done using `query_criteria`, where we specify:
\n", "\n", - "In this example, we align three F606W full-frame 339 second ACS/WFC images of the globular cluster NGC 104 from ACS/CAL Program [10737](http://www.stsci.edu/cgi-bin/get-proposal-info?id=10737&observatory=HST). These observations were acquired over a 3-month period in separate visits and at different orientations. We will use calibrated, CTE-corrected `*_flc.fits` files, which are available for both WFC3/UVIS and ACS/WFC detectors. \n", + "    --> obs_id, proposal_id, and filters \n", + "\n", + "MAST data products may be downloaded by using `download_products`, where we specify:
\n", + "\n", + "    --> products = calibrated (FLT, FLC) or drizzled (DRZ, DRC) files\n", + "\n", + "    --> type = standard products (CALxxx) or advanced products (HAP-SVM)\n", + "____\n", + "\n", + "\n", + "In this example, we align three F606W full-frame 339 second ACS/WFC images of the globular cluster NGC 104 from ACS/CAL Program [10737](http://www.stsci.edu/cgi-bin/get-proposal-info?id=10737&observatory=HST). These observations were acquired over a 3-month period in separate visits and at different orientations. We will use calibrated, CTE-corrected `*_flc.fits` files. \n", "\n", " j9irw3fwq_flc.fits\n", " j9irw4b1q_flc.fits\n", " j9irw5kaq_flc.fits\n", " \n", - "First, we use `astroquery` to look for the desired datasets and retrieve them from MAST. \n", - "\n", - "You may query on a large number of parameters, but to obtain these specific datasets we will only need to pass in a few: `obstype = 'all'` to include both calibration and GO programs in the search, `obs_collection = 'HST'` to exclude Hubble Legacy Archive images, and finally `obs_id` to search for a list of the dataset IDs. We finally select only 'FLC' data products using the parameter `productSubGroupDescription` while downloading the files. Because the second 'FLC' file for visit W4 is part of a dithered association, the `obs_id` for that assocation will need to be used instead of the individual filename. \n", - "\n", - "**Note**: The next cell may take awhile to complete or may need to be run more than once, as errors may occasionally arise from a failed server connection." + "
Depending on your connection speed, this cell may take a few minutes to execute.
" ] }, { @@ -79,27 +161,58 @@ "metadata": {}, "outputs": [], "source": [ - "# Query the data\n", - "obsTable = Observations.query_criteria(obstype='all', \n", - " obs_collection='HST',\n", - " obs_id=['j9irw3fwq', 'j9irw4040', 'j9irw5kaq'])\n", + "# Because the second 'FLC' file in visit W4 is part of a dithered association,\n", + "# the 'obs_id' for that assocation will need to be used instead of the individual filename.\n", + "\n", + "obs_ids = [\"j9irw3fwq\", \"j9irw4040\", \"j9irw5kaq\"]\n", + "props = [\"10737\"]\n", + "filts = [\"F606W\"]\n", "\n", - "# Download the files\n", + "obsTable = Observations.query_criteria(obs_id=obs_ids, proposal_id=props, filters=filts)\n", "products = Observations.get_product_list(obsTable)\n", - "Observations.download_products(products,download_dir='',\n", - " mrp_only=False,\n", - " productSubGroupDescription='FLC')\n", - "\n", - "# Move to working directory (not necessary, but outputs will all be in the same place if this is done)\n", - "input_flcs = glob.glob(os.path.join('mastDownload', 'HST', '*', '*flc.fits'))\n", - "print(input_flcs)\n", - "for flc in input_flcs:\n", - " shutil.copy(flc, os.path.basename(flc))\n", - "shutil.rmtree('mastDownload') #remove mast download dir now that we've moved the files\n", - "\n", - "#Remove the second file in the W4 association to simplify this notebook example and define the input file list.\n", - "os.remove('j9irw4b3q_flc.fits')\n", - "input_flcs = glob.glob('*flc.fits')\n", + "\n", + "data_prod = [\"FLC\"] # ['FLC','FLT','DRC','DRZ']\n", + "data_type = [\"CALACS\"] # ['CALACS','CALWF3','CALWP2','HAP-SVM']\n", + "\n", + "Observations.download_products(\n", + " products, productSubGroupDescription=data_prod, project=data_type, cache=True\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Move the files to the local working directory." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for flc in glob.glob(\"./mastDownload/HST/*/*flc.fits\"):\n", + " flc_name = os.path.split(flc)[-1]\n", + " os.rename(flc, flc_name)\n", + "shutil.rmtree(\"mastDownload/\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Remove the second FLC file in the W4 visit to simplify the example." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "os.remove(\"j9irw4b3q_flc.fits\")\n", + "input_flcs = glob.glob(\"*flc.fits\")\n", "print(input_flcs)" ] }, @@ -107,9 +220,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The files were downloaded into the subdirectory 'mastDownload' within the current working directory and then moved to the current working directory. \n", + "### 1.1 Check image header data \n", "\n", - "Let's take a look at the file headers to show the orientation differences between frames, as reflected by the 'PA_V3' keyword:" + "Here we will look at important keywords in the image headers" ] }, { @@ -118,35 +231,151 @@ "metadata": {}, "outputs": [], "source": [ - "for f in glob.glob('*flc.fits'):\n", - " print('filename:', fits.getval(f, 'rootname'), 'PA_V3:', fits.getval(f, 'PA_V3'))" + "paths = sorted(glob.glob(\"*flc.fits\"))\n", + "data = []\n", + "\n", + "ext_0_kws = [\n", + " \"rootname\",\n", + " \"ASN_ID\",\n", + " \"TARGNAME\",\n", + " \"detector\",\n", + " \"filter1\",\n", + " \"exptime\",\n", + " \"ra_targ\",\n", + " \"dec_targ\",\n", + " \"date-obs\",\n", + " \"postarg1\",\n", + " \"postarg2\",\n", + "]\n", + "ext_1_kws = [\"orientat\"]\n", + "\n", + "for path in paths:\n", + " path_data = []\n", + " for keyword in ext_0_kws:\n", + " path_data.append(fits.getval(path, keyword, ext=0))\n", + " for keyword in ext_1_kws:\n", + " path_data.append(fits.getval(path, keyword, ext=1))\n", + " data.append(path_data)\n", + "\n", + "keywords = ext_0_kws + ext_1_kws\n", + "table = Table(\n", + " np.array(data),\n", + " names=keywords,\n", + " dtype=[\n", + " \"str\",\n", + " \"str\",\n", + " \"str\",\n", + " \"str\",\n", + " \"str\",\n", + " \"f8\",\n", + " \"f8\",\n", + " \"f8\",\n", + " \"str\",\n", + " \"f8\",\n", + " \"f8\",\n", + " \"f8\",\n", + " ],\n", + ")\n", + "\n", + "table[\"exptime\"].format = \"7.1f\"\n", + "table[\"ra_targ\"].format = table[\"dec_targ\"].format = \"7.4f\"\n", + "table[\"postarg1\"].format = table[\"postarg2\"].format = \"7.3f\"\n", + "table[\"orientat\"].format = \"7.2f\"\n", + "table" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## 2. Align the Images with TweakReg\n", + "### 1.2 Inspect the Alignment \n", + "\n", + "Check the active WCS solution in the image header. If the image is aligned to a catalog, list the number of matches and the fit RMS in mas.
\n", + "Convert the fit RMS values to pixels for comparison with the alignment results performed later in this notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ext_0_kws = [\"DETECTOR\"]\n", + "ext_1_kws = [\"WCSNAME\", \"NMATCHES\", \"RMS_RA\", \"RMS_DEC\"]\n", "\n", - "In order to correctly align and combine multiple images into a distortion-free final drizzled product, `AstroDrizzle` relies on the World Coordinate System (WCS) information stored in the header of each image. Because the WCS information for an image is tied to the positions of guide stars used for the observation, each image will have small pointing errors offsets due to the uncertainty in guide star positions. The `TweakReg` software can align a set of images to better than subpixel accuracy (<0.1 pixel). All input images are aligned to one chosen reference image's WCS. The reference image is typically selected as the first image in the input list, but can be chosen explicitly via the `refimage` parameter, if desired.\n", + "det_scale = {\"IR\": 0.1283, \"UVIS\": 0.0396, \"WFC\": 0.05} # plate scale (arcsec/pixel)\n", "\n", - "The `TweakReg` algorithm consists of a few steps:\n", + "format_dict = {}\n", + "col_dict = defaultdict(list)\n", "\n", - " 1. Makes a catalog of source positions for each input (and reference) image using a source finding algorithm similar to DAOFind. \n", - " 2. Converts the pixel position catalogs to sky position catalogs.\n", - " 3. Finds common source positions between reference and input image catalogs.\n", - " 4. Calculates the shifts, rotation, and scale needed to align sky positions of sources in the input images.\n", - " 5. (Optionally) Updates the input image headers with the newly calculated WCS information." + "for f in sorted(glob.glob(\"*fl?.fits\")):\n", + " col_dict[\"filename\"].append(f)\n", + " hdr0 = fits.getheader(f, 0)\n", + " hdr1 = fits.getheader(f, 1)\n", + "\n", + " for kw in ext_0_kws: # extension 0 keywords\n", + " col_dict[kw].append(hdr0[kw])\n", + " for kw in ext_1_kws: # extension 1 keywords\n", + " if \"RMS\" in kw:\n", + " val = np.around(hdr1[kw], 1)\n", + " else:\n", + " val = hdr1[kw]\n", + " col_dict[kw].append(val)\n", + "\n", + " for kw in [\"RMS_RA\", \"RMS_DEC\"]:\n", + " val = np.round(\n", + " hdr1[kw] / 1000.0 / det_scale[hdr0[\"DETECTOR\"]], 2\n", + " ) # convert RMS from mas to pixels\n", + " col_dict[f\"{kw}_pix\"].append(val)\n", + "\n", + "wcstable = Table(col_dict)\n", + "wcstable" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we see that the FLC images from each separate visit have each been fit to the reference catalog GAIA eDR3. \n", + "\n", + "For visits w3 and w5, the WCSNAME 'FIT_IMG_catalog' implies that each FLC was aligned separately, whereas for visit w4, the drizzle combined image 'j9irw4040_drc.fits' was aligned to GAIA and the solution was propogated back to the FLC exposures. (Note that the number of matches and fit rms values are different for each FLC frame.)\n", + "\n", + "For more details on absolute astrometry in HST images, see [Section 4.5 in the DrizzlePac Handbook](https://hst-docs.stsci.edu/drizzpac/chapter-4-astrometric-information-in-the-header/4-5-absolute-astrometry).\n", + "\n", + "---" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Align with `TweakReg` \n", + "[Table of Contents](#toc)\n", + "\n", + "In order to combine multiple images into a distortion-free final drizzled product, `AstroDrizzle` relies on the World Coordinate System (WCS) information stored in the header of each image. Because the WCS information for an image is tied to the positions of guide stars used for the observation, each visit will have small pointing errors offsets due to the uncertainty in guide star positions. \n", + "\n", + "The `TweakReg` software can align a set of HST images to one another with subpixel accuracy (<0.1 pixel). All input images are aligned to one chosen reference image's WCS. The reference image is typically selected as the first image in the input list, but can be chosen explicitly via the `refimage` parameter, if desired.\n", + "\n", + "Note that since these images are already aligned to GAIA with ~400 matches and fit RMS values ~0.2 pixels. Therefore, we simply want to check the relative alignment of the FLC frame in the individual visits W3, W4 and W5, building on the eDR3 WCS solution and then testing for any residual offsets between frames. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ + "The `TweakReg` algorithm consists of a few steps:\n", + "\n", + " 1. Makes a catalog of source positions for each input (and reference) image using a source finding algorithm similar to DAOFind. \n", + " 2. Converts the pixel position catalogs to sky position catalogs.\n", + " 3. Finds common source positions between reference and input image catalogs.\n", + " 4. Calculates the shifts, rotation, and scale needed to align sky positions of sources in the input images.\n", + " 5. (Optionally) Updates the input image headers with the newly calculated WCS information.\n", + " \n", + " \n", "`TweakReg` has numerous parameters that are listed in the documentation. For this first pass, we will run `TweakReg` with default parameters and show how it fails for this dataset. In subsequent attempts, we will change critical parameters and show how this improves the alignment. The 5 parameters listed below are commented out when left at their default values and uncommented when set to non-defaults. \n", "\n", - " 'conv_width' : 2x FWHM (~3.5 pix for point sources in ACS/WFC or WFC3/UVIS)\n", - " 'threshold' : Signal-to-noise above background. Start large and reduce until number of objects is acceptable.\n", + " 'conv_width' : 2x FWHM (~3.5 pix for point sources in ACS/WFC or WFC3/UVIS, ~2.5 pix for WFC3/IR)\n", + " 'threshold' : Signal-to-noise above bkg. Start large and reduce until number of objects is acceptable.\n", " 'peakmax' : Source brightness cutoff, set to avoid saturated sources. \n", " 'searchrad' : Search radius for finding common sources between images.\n", " 'fitgeometry': Geometry used to fit offsets, rotations and/or scale changes from the matched object lists. \n", @@ -160,8 +389,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### 2a. First test: Use 'default' parameters. \n", - "### NOTE: DONT ACTUALLY RUN THE CELLS IN THIS EXAMPLE, THEY USE TOO MUCH MEMORY!\n", + "### 2.1 First test: Use 'default' parameters. \n", + "\n", + "
NOTE: Do not actually run the cells in this example, as they will use too much memory!
\n", "\n", "The most basic, out-of-the-box call to TweakReg looks like this. You `TweakReg` on your list of input files." ] @@ -188,9 +418,8 @@ " > 'outshifts' is the name of the output file created when `shiftfile` is set to True.\n", " \n", " > 'updatehdr' is False by default, and specifies whether or not to update the headers of each input image directly with the shifts that were determined. This will allow the input images to be combined by AstroDrizzle without having to provide the shiftfile as well.\n", - " \n", " \n", - "There are many more parameters that can be set, but the example below shows how to run `TweakReg` and override some of these defaults.\n" + "There are many more parameters that can be set, but the example below shows how to run `TweakReg` and override some of these defaults." ] }, { @@ -199,27 +428,52 @@ "metadata": {}, "outputs": [], "source": [ - "# running tweakreg with some parameters changed from their defaults.\n", + "# Running tweakreg with some parameters changed from their defaults.\n", + "#\n", "# tweakreg.TweakReg(input_flcs,\n", "# interactive=True,\n", - "# shiftfile=True, \n", - "# outshifts='shift_default.txt', \n", - "# updatehdr=True)" + "# shiftfile=True,\n", + "# outshifts='shift_default.txt',\n", + "# updatehdr=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### 2b. Second test: Setting parameters for source finding.\n", + "### 2.2 Second test: Set parameters for source finding \n", "\n", "One of the crucial steps for aligning images is creating a catalog of common sources in each image to do the alignment. `TweakReg` uses a set of default parameters to do this source finding, but you will likely want to override these and tailor them to your dataset. As described above, using the defaults on this dataset of a cluster will simply produce too many sources and use too much memory for the average user's machine. In this example, we will adjust some of these to parameters.\n", "\n", "In the above example, none of the source finding parameters were adjusted. `TweakReg`, by default, looks for sources which are 4-sigma above the background, but this `threshold` value is too low, since the catalogs contain ~60,000 objects per chip in each frame, and ~6500 matched sources between each exposure. Setting the `threshold` to a very low value generally does *not* translate to a better solution, since all sources are weighted equally when computing fits to match catalogs. This is especially relevant for ACS/WFC and WFC3/UVIS data where CTE tails can shift the centroid position slightly along the readout direction for faint sources and potentially bias the fit.\n", "\n", - "In this test, the `threshold` value is adjusted to a larger value, and the `peakmax` value is set to 70,000 electrons to exclude full-well saturated objects for gain=2. `TweakReg` now finds a more reasonable number of sources per ACS chip (~500).\n", + "In this test, the `threshold` value is adjusted to a larger value, and the `peakmax` value is set to 70,000 electrons to exclude full-well saturated objects for gain=2. `TweakReg` now finds a more reasonable number of sources per ACS chip (~700).\n", + "\n", + "Note that there are other source finding parameters that can be adjusted and that you may need to adjust. For example, you may need to also adjust `searchrad` if your images have large offsets." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "tweakreg.TweakReg(\n", + " input_flcs,\n", + " threshold=4000,\n", + " peakmax=70000,\n", + " configobj=None,\n", + " interactive=False,\n", + " shiftfile=True,\n", + " outshifts=\"shift_rscale.txt\",\n", + " ylimit=0.2,\n", + " fitgeometry=\"rscale\",\n", + " updatehdr=False,\n", + ")\n", "\n", - "Note that there are other source finding paratmeters that can be adjusted and that you may need to adjust. For example, you may need to also adjust `searchrad` up if your images have large offsets." + "clear_output()" ] }, { @@ -228,26 +482,25 @@ "metadata": {}, "outputs": [], "source": [ - "tweakreg.TweakReg(input_flcs, \n", - " threshold=4000, \n", - " peakmax=70000,\n", - " configobj = None, \n", - " interactive=False,\n", - " shiftfile=True,\n", - " outshifts='shift_thresh.txt', \n", - " updatehdr=False)" + "# If the alignment is unsuccessful, stop the notebook\n", + "with open(\"shift_rscale.txt\", \"r\") as shift:\n", + " for line_number, line in enumerate(shift, start=1):\n", + " if \"nan\" in line:\n", + " raise ValueError(\"nan found in line {} in shift file\".format(line_number))\n", + " else:\n", + " continue" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## 2c. Inspecting the quality of the fit\n", + "### 2.3 Inspect the quality of the fit \n", "\n", "In your work, you will likely have to play around with TweakReg parameters to find the optimal set for your data. This requires inspecting the quality of the fit after each run. If you scroll to the bottom of the last cell run, you can see some diagnostic plots there including a residual and vector plot that show the quality of fit. Here, we will inspect some of the outputs that are saved.\n", "\n", "\n", - "The shift file below shows that the xrms and yrms of the computed fit are decent, around ~0.04 pixels." + "The shift file below shows that the xrms and yrms of the computed fit are excellent, around ~0.04 pixels." ] }, { @@ -256,11 +509,13 @@ "metadata": {}, "outputs": [], "source": [ - "shift_tab = Table.read('shift_thresh.txt',\n", - " format='ascii.no_header',\n", - " names=['file','dx','dy','rot','scale','xrms','yrms'])\n", + "shift_tab = Table.read(\n", + " \"shift_rscale.txt\",\n", + " format=\"ascii.no_header\",\n", + " names=[\"file\", \"dx\", \"dy\", \"rot\", \"scale\", \"xrms\", \"yrms\"],\n", + ")\n", "\n", - "formats = ['.2f', '.2f', '.3f', '.5f', '.2f', '.2f']\n", + "formats = [\".2f\", \".2f\", \".3f\", \".5f\", \".2f\", \".2f\"]\n", "for i, col in enumerate(shift_tab.colnames[1:]):\n", " shift_tab[col].format = formats[i]\n", "shift_tab" @@ -279,7 +534,21 @@ "metadata": {}, "outputs": [], "source": [ - "Image(filename='residuals_j9irw5kaq_flc.png', width=500, height=600) " + "# Give the 'fit residual plots' a unique name for comparison with other tests.\n", + "rootname_A = \"j9irw4b1q\"\n", + "rootname_B = \"j9irw5kaq\"\n", + "\n", + "# make a subdirectory for storing images\n", + "try:\n", + " os.mkdir(\"images\")\n", + "except FileExistsError:\n", + " pass\n", + "\n", + "shutil.copy(f\"residuals_{rootname_A}_flc.png\", f\"images/residuals_{rootname_A}_flc_rscale.png\")\n", + "shutil.copy(f\"residuals_{rootname_B}_flc.png\", f\"images/residuals_{rootname_B}_flc_rscale.png\")\n", + "\n", + "shutil.copy(f\"vector_{rootname_A}_flc.png\", f\"images/vector_{rootname_A}_flc_rscale.png\")\n", + "shutil.copy(f\"vector_{rootname_B}_flc.png\", f\"images/vector_{rootname_B}_flc_rscale.png\")" ] }, { @@ -288,33 +557,67 @@ "metadata": {}, "outputs": [], "source": [ - "Image(filename='vector_j9irw4b1q_flc.png', width=500, height=600) " + "# read images\n", + "img_A = mpimg.imread(f\"images/residuals_{rootname_A}_flc_rscale.png\")\n", + "img_B = mpimg.imread(f\"images/residuals_{rootname_B}_flc_rscale.png\")\n", + "\n", + "# display images\n", + "fig, axs = plt.subplots(1, 2, figsize=(10, 10), dpi=200)\n", + "axs[0].imshow(img_A)\n", + "axs[1].imshow(img_B)\n", + "axs[0].axis(\"off\")\n", + "axs[1].axis(\"off\")\n", + "fig.tight_layout()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# read images\n", + "img_A = mpimg.imread(f\"images/vector_{rootname_A}_flc_rscale.png\")\n", + "img_B = mpimg.imread(f\"images/vector_{rootname_B}_flc_rscale.png\")\n", + "\n", + "# display images\n", + "fig, axs = plt.subplots(1, 2, figsize=(10, 10), dpi=100)\n", + "axs[0].imshow(img_A)\n", + "axs[1].imshow(img_B)\n", + "axs[0].axis(\"off\")\n", + "axs[1].axis(\"off\")\n", + "fig.tight_layout()\n", + "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The residual plots show a decent RMS of fit around 0.04 pixels, but there is still a systematic skew in the fit residuals, which is believed to be caused by an uncorrected 'skew' in the ACS distortion solution. \n", + "The residual plots show a fit RMS of ~0.04 pixels, but there is still a systematic skew in the fit residuals, which is believed to be caused by an uncorrected 'skew' in the ACS distortion solution. \n", + "\n", + "As of creation of this notebook in 2024, the following calibration reference files were used:\n", "\n", - "As of creation of this notebook in 2018, the following calibration reference files were used:\n", + "- IDCTAB = 'jref$4bb1536oj_idc.fits' / Image Distortion Correction Table\n", "\n", - " IDCTAB = 'jref$0461802ej_idc.fits' / Image Distortion Correction Table\n", - " D2IMFILE = 'jref$02c1450oj_d2i.fits' / Column Correction Reference File\n", - " NPOLFILE = 'jref$02c14514j_npl.fits' / Non-polynomial Offsets Reference File (F606W)\n", + "- D2IMFILE = 'jref$4bb15371j_d2i.fits' / Column Correction Reference File \n", + "\n", + "- NPOLFILE = 'jref$4bb1536ej_npl.fits' / Non-polynomial Offsets Reference File (F606W) \n", " \n", - "New distortion solutions, based on the latest Gaia DR2, are in the process of being derived by the ACS team, so observations retrieved from MAST at a later date may or may not show these systematic skew residuals. Fortunately, these can be corrected by allowing for a higher order `fitgeometry` as shown in the next test.\n", + "New distortion solutions are in the process of being derived by the ACS team, so observations retrieved from MAST at a later date may or may not show these systematic skew residuals. Fortunately, these can be corrected by allowing for a higher order `fitgeometry` as shown in the next test.\n", "\n", - "More information about the distortion reference files used by `AstroDrizzle` may be found [here](http://documents.stsci.edu/hst/HST_overview/documents/DrizzlePac/ch33.html)." + "More information about the distortion reference files used by `AstroDrizzle` may be found in the DrizzlePac Handbook [Section 4.3](https://hst-docs.stsci.edu/drizzpac/chapter-4-astrometric-information-in-the-header/4-3-distortion-information-in-pipeline-calibrated-images)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### 2d. Fourth test: Adjust the `fitgeometry`\n", + "### 2.4 Fourth test: Adjust the `fitgeometry` \n", "\n", - "In this test, the `fitgeometry` parameter is changed from the default 'rscale' to 'general' in order to allows for an xy/shift and an independent scale and rotation for each axis." + "In this test, the `fitgeometry` parameter is changed from the default `rscale` (shift, rotation, scale) to `general` (shift, xrot, yrot, xscale, yscale) in order to correct for residual skew terms." ] }, { @@ -323,16 +626,19 @@ "metadata": {}, "outputs": [], "source": [ - "tweakreg.TweakReg(input_flcs, \n", - " threshold=4000, \n", - " searchrad=4.0,\n", + "tweakreg.TweakReg(\n", + " input_flcs,\n", + " threshold=4000,\n", " peakmax=70000,\n", - " fitgeometry='general',\n", - " configobj = None, \n", + " configobj=None,\n", " interactive=False,\n", - " shiftfile=True, \n", - " outshifts='shift_general.txt', \n", - " updatehdr=False)" + " shiftfile=True,\n", + " outshifts=\"shift_general.txt\",\n", + " ylimit=0.2,\n", + " fitgeometry=\"general\",\n", + " updatehdr=False,\n", + ")\n", + "clear_output()" ] }, { @@ -341,26 +647,57 @@ "metadata": {}, "outputs": [], "source": [ - "shift_tab=Table.read('shift_general.txt',\n", - " format='ascii.no_header',\n", - " names=['file','dx','dy','rot','scale','xrms','yrms'])\n", + "shift_tab = Table.read(\n", + " \"shift_general.txt\",\n", + " format=\"ascii.no_header\",\n", + " names=[\"file\", \"dx\", \"dy\", \"rot\", \"scale\", \"xrms\", \"yrms\"],\n", + ")\n", "\n", - "formats = ['.2f', '.2f', '.3f', '.5f', '.2f', '.2f']\n", + "formats = [\".2f\", \".2f\", \".3f\", \".5f\", \".2f\", \".2f\"]\n", "for i, col in enumerate(shift_tab.colnames[1:]):\n", " shift_tab[col].format = formats[i]\n", "shift_tab" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Print the number of matches for images 2 and 3 with respect to image1\n", + "rootname2 = \"j9irw4b1q\"\n", + "rootname3 = \"j9irw5kaq\"\n", + "match_tab2 = ascii.read(\n", + " rootname2 + \"_flc_catalog_fit.match\"\n", + ") # load match file in astropy table\n", + "match_tab3 = ascii.read(\n", + " rootname3 + \"_flc_catalog_fit.match\"\n", + ") # load match file in astropy table\n", + "print(rootname2, \" Number of Gaia Matches=\", {len(match_tab2)})\n", + "print(rootname3, \" Number of Gaia Matches=\", {len(match_tab3)})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Compare with the MAST solutions\n", + "wcstable" + ] + }, { "cell_type": "markdown", "metadata": {}, "source": [ - "While the shift file looks nearly identical to the prior run, this is because the axis-dependent rotation and scale values are averaged together in the shiftfile. To determine the actual values, it is necessary to inspect the logfile.\n", + "While the shift file looks nearly identical to the prior run, this is because the axis-dependent rotation and scale values are averaged together in the shiftfile. To determine the actual values, it is necessary to inspect the 'tweakreg.log' logfile.\n", "\n", - " j9irw5kaq SCALE_X: 1.000008792 SCALE_Y: 1.00000757 ROT_X: 0.002818152 ROT_Y: 359.998708 SKEW: 359.9958898\n", - " j9irw4b1q SCALE_X: 0.999986877 SCALE_Y: 1.00003072 ROT_X: 0.004877188 ROT_Y: 0.00405990 SKEW: -0.0008173 \n", + " j9irw4b1q SCALE_X: 0.9999720117 SCALE_Y: 1.000014811 ROT_X: 359.9999 ROT_Y: 0.0014 SKEW: -359.9985\n", + " j9irw5kaq SCALE_X: 0.9999851345 SCALE_Y: 1.000030317 ROT_X: 0.0015 ROT_Y: 359.9989 SKEW: 359.9974\n", "\n", - "The 'general' fit reduces the fit rms to ~0.03 pixels from the prior value of ~0.04 pixels and removes systematics from the residuals. " + "The 'general' fit reduces the fit rms to ~0.03 pixels from the prior value of ~0.05 pixels and removes systematics from the residuals in the plots below. " ] }, { @@ -369,17 +706,39 @@ "metadata": {}, "outputs": [], "source": [ - "#Give the 'fit residual plots' a unique name for comparison with other tests.\n", + "# Give the 'fit residual plots' a unique name for comparison with other tests.\n", + "shutil.copy(f\"residuals_{rootname_A}_flc.png\", f\"images/residuals_{rootname_A}_flc_general.png\")\n", + "shutil.copy(f\"residuals_{rootname_B}_flc.png\", f\"images/residuals_{rootname_B}_flc_general.png\")\n", "\n", - "os.rename('residuals_j9irw4b1q_flc.png', 'residuals_j9irw4b1q_flc_general.png')\n", - "Image(filename='residuals_j9irw4b1q_flc_general.png', width=500, height=600) " + "shutil.copy(f\"vector_{rootname_A}_flc.png\", f\"images/vector_{rootname_A}_flc_general.png\")\n", + "shutil.copy(f\"vector_{rootname_B}_flc.png\", f\"images/vector_{rootname_B}_flc_general.png\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The plot above shows residuals in X and Y plotted as functions of X and Y. Each point represents a source that was used for the alignment. The residuals should look fairly random - if any correlation is seen, this is an indicator of a poor alignment solution. In some cases you will notice a wavy pattern in the residuals in X and Y. This is often seen for UVIS images and is a result of lithographic patterns of the detector that are not fully corrected for in the distortion solutions. The RMS in X and Y are also printed. For a good alignment, we are looking for an RMS on the order of 0.1 pixel or less." + "The plots below shows residuals in X and Y plotted as functions of X and Y. Each point represents a source that was used for the alignment. The residuals should look fairly random - if any correlation is seen, this is an indicator of a poor alignment solution. In some cases you will notice a wavy pattern in the residuals in X and Y. This is often seen for UVIS images and is a result of lithographic patterns of the detector that are not fully corrected for in the distortion solutions. The RMS in X and Y are also printed. For a good alignment, we are looking for an RMS on the order of 0.1 pixel or less." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# read images\n", + "img_A = mpimg.imread(f\"images/residuals_{rootname_A}_flc_general.png\")\n", + "img_B = mpimg.imread(f\"images/residuals_{rootname_B}_flc_general.png\")\n", + "\n", + "# display images\n", + "fig, axs = plt.subplots(1, 2, figsize=(20, 10), dpi=200)\n", + "axs[0].imshow(img_A)\n", + "axs[1].imshow(img_B)\n", + "axs[0].axis(\"off\")\n", + "axs[1].axis(\"off\")\n", + "fig.tight_layout()\n", + "plt.show()" ] }, { @@ -388,10 +747,18 @@ "metadata": {}, "outputs": [], "source": [ - "#Give the 'vector residual plots' a unique name for comparison with other tests.\n", + "# read images\n", + "img_A = mpimg.imread(f\"images/vector_{rootname_A}_flc_general.png\")\n", + "img_B = mpimg.imread(f\"images/vector_{rootname_B}_flc_general.png\")\n", "\n", - "os.rename('vector_j9irw4b1q_flc.png', 'vector_j9irw4b1q_flc_general.png')\n", - "Image(filename='vector_j9irw4b1q_flc_general.png', width=500, height=600) " + "# display images\n", + "fig, axs = plt.subplots(1, 2, figsize=(10, 10), dpi=100)\n", + "axs[0].imshow(img_A)\n", + "axs[1].imshow(img_B)\n", + "axs[0].axis(\"off\")\n", + "axs[1].axis(\"off\")\n", + "fig.tight_layout()\n", + "plt.show()" ] }, { @@ -412,9 +779,11 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### 2e. Update header once optimal parameters are found\n", + "### 2.5 Overplot Matched Sources on the Image \n", + "\n", + "As a final quality check before updating the image header WCS, we plot the sources that were matched between FLC images in visits W3 and W4.\n", "\n", - "Once you've decided on the optimal set of parameters for your alignment, it is safe to update the header of the input data." + "The cell below shows how to read in the `*_fit.match` file as an `astropy` table. Unfortunately, it doesn't automatically name columns so you'll have to look at the header to grab the columns." ] }, { @@ -423,27 +792,50 @@ "metadata": {}, "outputs": [], "source": [ - "# Final run with ideal parameters, updatehdr = True\n", + "rootname = \"j9irw4b1q\"\n", "\n", - "tweakreg.TweakReg(input_flcs, \n", - " threshold=4000, \n", - " peakmax=70000,\n", - " fitgeometry='general',\n", - " configobj = None, \n", - " interactive=False,\n", - " shiftfile=False, \n", - " updatehdr=True)" + "plt.figure(figsize=(7, 7), dpi=140)\n", + "chip1_data = fits.open(rootname + \"_flc.fits\")[\"SCI\", 2].data\n", + "chip2_data = fits.open(rootname + \"_flc.fits\")[\"SCI\", 1].data\n", + "fullsci = np.concatenate([chip2_data, chip1_data])\n", + "zscale = ZScaleInterval()\n", + "z1, z2 = zscale.get_limits(fullsci)\n", + "plt.imshow(fullsci, cmap=\"Greys\", origin=\"lower\", vmin=z1, vmax=z2)\n", + "\n", + "match_tab = ascii.read(\n", + " rootname + \"_flc_catalog_fit.match\"\n", + ") # load match file in astropy table\n", + "match_tab_chip1 = match_tab[\n", + " match_tab[\"col15\"] == 2\n", + "] # filter table for sources on chip 1 (on ext 4)\n", + "match_tab_chip2 = match_tab[\n", + " match_tab[\"col15\"] == 1\n", + "] # filter table for sources on chip 1 (on ext 4)\n", + "x_cord1, y_cord1 = match_tab_chip1[\"col11\"], match_tab_chip1[\"col12\"]\n", + "x_cord2, y_cord2 = match_tab_chip2[\"col11\"], match_tab_chip2[\"col12\"]\n", + "\n", + "plt.scatter(\n", + " x_cord1,\n", + " y_cord1 + 2051,\n", + " s=50,\n", + " edgecolor=\"r\",\n", + " facecolor=\"None\",\n", + " label=\"Matched Sources\",\n", + ")\n", + "plt.scatter(x_cord2, y_cord2, s=50, edgecolor=\"r\", facecolor=\"None\")\n", + "plt.title(\n", + " f\"Matched sources W4 to W3: N = {len(match_tab)}\", fontsize=14\n", + ")\n", + "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### 2f. Overplot Matched Sources on the Image\n", - "\n", - "Let's plot the sources that were matched between all images on the bottom chip. Confusingly, this is referred to as 'chip 2' or SCI,1 or extension 1. \n", + "### 2.6 Update header once optimal parameters are found \n", "\n", - "The cell below shows how to read in the `*_fit.match` file as an `astropy` table. Unfortunatley, it doesn't automatically name columns so you'll have to look at the header to grab the columns." + "Once you've decided on the optimal set of parameters for your alignment, it is safe to update the header of the input data. Note that because the first image is already aligned to Gaia, this does not need to be done again and we focus only on improving the relative alignment." ] }, { @@ -452,27 +844,27 @@ "metadata": {}, "outputs": [], "source": [ - "plt.figure(figsize = (20, 10))\n", - "chip1_data = fits.open('j9irw4b1q_flc.fits')['SCI', 1].data\n", - "zscale = ZScaleInterval()\n", - "z1, z2 = zscale.get_limits(chip1_data)\n", - "plt.imshow(chip1_data, cmap='Greys',origin='lower', vmin=z1, vmax=z2)\n", - "\n", - "match_tab = ascii.read('j9irw4b1q_flc_catalog_fit.match') #load match file in astropy table\n", - "match_tab_chip2 = match_tab[match_tab['col15'] == 1] #filter table for sources on chip 2 (on ext 1)\n", - "x_cord, y_cord = match_tab_chip2['col11'], match_tab_chip2['col12']\n", + "tweakreg.TweakReg(\n", + " input_flcs,\n", + " threshold=4000,\n", + " peakmax=70000,\n", + " configobj=None,\n", + " interactive=False,\n", + " shiftfile=False,\n", + " ylimit=0.2,\n", + " fitgeometry=\"general\",\n", + " updatehdr=True, # update the header with the new solution\n", + ") \n", "\n", - "plt.scatter(x_cord, y_cord, s=30, edgecolor='r', facecolor='None',label='Matched Sources, Chip 2')\n", - "plt.ylim(0,2051)\n", - "plt.xlim(0,4096)\n", - "plt.legend(loc='best', fontsize=20)" + "clear_output()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## 3. Combine the Images using AstroDrizzle\n", + "## 3. Combine the Images using `AstroDrizzle` \n", + "[Table of Contents](#toc)\n", "\n", "Now that the images are aligned to a common WCS, they are ready for combination with `AstroDrizzle`. All of the input exposures taken in a single filter will contribute to a single drizzled output file.\n", "\n", @@ -483,8 +875,7 @@ " 3. Each image is individually drizzled, with geometric distortion corrections, to a common reference frame.\n", " 4. The distortion-free drizzled images are combined to create a median image.\n", " 5. The median image is blotted, or reverse-drizzled, back to the frame of each input image.\n", - " 6. By comparing each input image with its counterpart blotted median image, the software locates bad pixels in\n", - " each of the original frames & creates bad pixel masks (typically cosmic rays and bad pixels in the detector)\n", + " 6. By comparing each input image with its counterpart blotted median image, the software locates bad pixels in each of the original frames & creates bad pixel masks (typically CR's and bad pixels in the detector)\n", " 7. In the final step, input images are drizzled together onto a single output image." ] }, @@ -494,7 +885,55 @@ "source": [ "Let's first run `AstroDrizzle` on the aligned input images in the next cell to create a combined image `f606w_combined_drc`, then go into further detail about the drizzle process. \n", "\n", - "Note that astrodrizzle supports the TEAL GUI interface for setting parameters, as well as loading in a custom configuration file (`.cfg`) files, but we will be using the command-line syntax interface to set the parameters in this example, where parameters are passed into the function directly. Any existing `.cfg` file will be overridden by setting `configobj = None` so that unless explicitly set, parameters will be reset to default." + "Note that astrodrizzle supports the loading of a custom configuration file (`.cfg`) files, but we will be using the command-line syntax interface to set the parameters in this example, where parameters are passed into the function directly. Any existing `.cfg` file will be overridden by setting `configobj = None` so that unless explicitly set, parameters will be reset to default." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, we get some recommended values for drizzling from the MDRIZTAB reference file. The parameters in this file are different for each detector and are based on the number of input frames. These are a good starting point for drizzling and may be adjusted accordingly." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# The following lines of code find and download the MDRIZTAB reference file.\n", + "mdz = fits.getval(input_flcs[0], 'MDRIZTAB', ext=0).split('$')[1]\n", + "print('Searching for the MDRIZTAB file:', mdz)\n", + "get_mdriztab = os.system('crds sync --hst --files ' + mdz + ' --output-dir '+os.environ['jref'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def get_vals_from_mdriztab(\n", + " input_flcs,\n", + " kw_list=[\n", + " \"driz_sep_bits\",\n", + " \"combine_type\",\n", + " \"driz_cr_snr\",\n", + " \"driz_cr_scale\",\n", + " \"final_bits\",\n", + " ],\n", + "):\n", + " \"\"\"Get only selected parameters from. the MDRIZTAB\"\"\"\n", + " mdriz_dict = getMdriztabPars(input_flcs)\n", + "\n", + " requested_params = {}\n", + "\n", + " print(\"Outputting the following parameters:\")\n", + " for k in kw_list:\n", + " requested_params[k] = mdriz_dict[k]\n", + " print(k, mdriz_dict[k])\n", + "\n", + " return requested_params" ] }, { @@ -503,15 +942,45 @@ "metadata": {}, "outputs": [], "source": [ - "astrodrizzle.AstroDrizzle(input_flcs,\n", - " output='f606w_combined',\n", + "selected_params = get_vals_from_mdriztab(input_flcs)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that the parameter `final_bits`= '336' is equivalent to `final_bits`= '256, 64, 16'. \n", + "\n", + "Since the ACS team now corrects for stable hot pixels (DQ flag=16) via the dark reference files, these pixels can be considered 'good'. Full well saturated pixels (DQ flag=256) and warm pixels (DQ flag=64) may also be treated as good. More details on the recommended drizzle parameters for ACS may be found in [ISR 2017-02](https://ui.adsabs.harvard.edu/abs/2017acs..rept....2H/abstract).\n", + "\n", + "Next we run `AstroDrizzle` to remove the sky, flag cosmic rays, and combine the image using the selected parameters. Note that these may be further refined by uncommenting the lines in the example shown below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# To override any of the above values, e.g.:\n", + "# selected_params[\"driz_sep_bits\"] = \"256, 64, 16\"\n", + "# selected_params[\"final_bits\"] = \"256, 64, 16\"\n", + "\n", + "astrodrizzle.AstroDrizzle(\n", + " input_flcs,\n", + " output=\"f606w_combined\",\n", " preserve=False,\n", - " driz_sep_bits='64,16',\n", + " build=True,\n", + " skymethod=\"match\",\n", + " driz_sep_bits=selected_params[\"driz_sep_bits\"],\n", " driz_cr_corr=True,\n", - " final_bits='64,16',\n", " clean=False,\n", + " final_bits=selected_params[\"final_bits\"],\n", " configobj=None,\n", - " build=True)" + ")\n", + "clear_output()" ] }, { @@ -529,7 +998,7 @@ " \n", " When `build = False`, these will be output as separate files (`_sci.fits`, `_wht.fits`, `_ctx.fits`).\n", "\n", - " When `final_wht_type = EXP` (default), the weight image it is effectively an exposure time map of each pixel in the final drizzled image. Other options are error 'ERR' and inverse variance 'IVM' weighting, as described [in more detail here](https://drizzlepac.readthedocs.io/en/latest/astrodrizzle.html).\n", + " When `final_wht_type = EXP` (default), the weight image it is effectively an exposure time map of each pixel in the final drizzled image. Other options are error 'ERR' and inverse variance 'IVM' weighting, as described [in more detail here](https://drizzlepac.readthedocs.io/en/latest/drizzlepac_api/astrodrizzle.html).\n", " \n", " The context image is a map showing which images contribute to the final drizzled stack. Each input image chip is identified by a bit in a 32-bit integer. For example, image1/chip1 = 2^0 = 1, image1/chip2 = 2^1 = 2. Each context image pixel is an additive combination of these bits, depending on which images contributed to the corresponding pixel in the drizzled image.\n", " \n", @@ -549,7 +1018,16 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Inspect the final drizzled image" + "### 3.1 Inspect the drizzled science image " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, we inspect the science and weight images and inspect their data quality. An imprint of sources in the weight image may indicate an error with the alignment and subsequent cosmic-ray rejection. \n", + "\n", + "For more details on inspecting the drizzled products after reprocessing, see [Section 7.3 in the DrizzlePac Handbook](https://hst-docs.stsci.edu/drizzpac/chapter-7-data-quality-checks-and-trouble-shooting-problems/7-3-inspecting-drizzled-products-after-user-reprocessing)" ] }, { @@ -558,18 +1036,21 @@ "metadata": {}, "outputs": [], "source": [ - "plt.figure(figsize = (10, 10))\n", - "drc_dat = fits.open('f606w_combined_drc.fits')['SCI', 1].data #final drizzled image in SCI,1 extension\n", + "plt.figure(figsize=(8, 8))\n", + "drc_dat = fits.open(\"f606w_combined_drc.fits\")[\n", + " \"SCI\", 1\n", + "].data # final drizzled image in SCI,1 extension\n", "z1, z2 = zscale.get_limits(drc_dat)\n", - "plt.imshow(drc_dat, origin='lower', vmin=z1, vmax=z2, cmap='Greys')\n", - "plt.title('F606W drizzled image', fontsize=30)" + "plt.imshow(drc_dat, origin=\"lower\", vmin=z1, vmax=z2, cmap=\"Greys\")\n", + "plt.title(\"F606W drizzled science image\", fontsize=20)\n", + "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Inspect the final weight image" + "### 3.2 Inspect the final weight image " ] }, { @@ -578,18 +1059,21 @@ "metadata": {}, "outputs": [], "source": [ - "plt.figure(figsize = (10, 10))\n", - "drc_dat = fits.open('f606w_combined_drc.fits')['WHT', 1].data #final drizzled image in WHT,1 extension\n", + "plt.figure(figsize=(8, 8))\n", + "drc_dat = fits.open(\"f606w_combined_drc.fits\")[\n", + " \"WHT\", 1\n", + "].data # final drizzled image in WHT,1 extension\n", "z1, z2 = zscale.get_limits(drc_dat)\n", - "plt.imshow(drc_dat, origin='lower', vmin=z1, vmax=z2, cmap='Greys')\n", - "plt.title('F606W weight image', fontsize=30)" + "plt.imshow(drc_dat, origin=\"lower\", vmin=z1, vmax=z2, cmap=\"Greys\")\n", + "plt.title(\"F606W drizzled weight image\", fontsize=20)\n", + "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Discussion of AstroDrizzle" + "### 3.3 Discussion of AstroDrizzle and Summary " ] }, { @@ -598,29 +1082,66 @@ "source": [ "In this example, we imported and ran the `AstroDrizzle` task on our input images, which were previously aligned with `TweakReg`. By setting `configobj` to None, we ensured that `AstroDrizzle` was not picking up any existing configuration files and parameters were restored to default values. We then set a select few parameters to non-default values:\n", " 1. `output = 'f606w_combined'` : Output file name root (which will be appended with various suffixes). Defaults to input file name.\n", - " 2. `driz_sep_bits = '64,16'` : Data quality flags in the `flt.fits` file, which were set during calibration, can be used as bit mask when drizzling. The user may specify which bit values should actually be considered \"good\" and included in image combination. In `astrodrizzle`, this parameter may be given as the sum of those DQ flags or as a comma-separated list, as shown in this example. In this example, 64 and 16 are set, so that both warm pixels and stable hot pixels, which are corrected by the dark, are treated as valid input pixels.\n", + " 2. `driz_sep_bits = '256, 64, 16'` : Data quality flags in the `flt.fits` file, which were set during calibration, can be used as bit mask when drizzling. The user may specify which bit values should actually be considered \"good\" and included in image combination. In `astrodrizzle`, this parameter may be given as the sum of those DQ flags or as a comma-separated list, as shown in this example. In this example, 64 and 16 are set, so that both warm pixels and stable hot pixels, which are corrected by the dark, are treated as valid input pixels. To avoid flagging saturated pixels as CR's, dq 256 is also set as a valid input pixel.\n", " 3. `driz_cr_corr = True` : When set to True, the task will create both a cosmic ray mask image (suffix `crmask.fits`) and a clean version of the original input images (suffix `crclean.fits`), where flagged pixels are replaced by pixels from the blotted median. It is strongly recommended that the quality of the cosmic ray masks be verified by blinking the original `flt.fits` input image with both the cosmic ray-cleaned image (`crclean.fits`) and the cosmic-ray mask (`crmask.fits`).\n", - " 4. `final_bits = '64,16'` : Similar to `driz_sep_bits`, but for the last step when all the input images are combined.\n", + " 4. `final_bits = '256, 64, 16'` : Similar to `driz_sep_bits`, but for the last step when all the input images are combined. To avoid filling the cores of saturated pixels with NaN, dq 256 is set as a valid input pixel.\n", " 5. `clean = False` : intermediate output files (e.g masks) will be kept. If True, only main outputs will be kept. Also see `in_memory` to control this behavior.\n", - " 6. `configobj = None` : ignore any TEAL inputs / config files and refresh all parameters to default values.\n", + " 6. `configobj = None` : ignore any custom config files and refresh all parameters to default values.\n", " \n", - "`AstroDrizzle` has a large number of parameters, which are [described here](https://drizzlepac.readthedocs.io/en/deployment/astrodrizzle.html). Running `AstroDrizzle` using default parameter values is not recommended, as these defaults may not provide optimal science products. Users should also inspect the quality of the sky subtraction and cosmic ray rejection. For dithered data, users may experiment with the output `final_scale` and `final_pixfrac` parameters in the `final_drizzle` step. For more details, see the notebook in this repository to 'Optimize Image Sampling'." + "`AstroDrizzle` has a large number of parameters, which are [described here](https://drizzlepac.readthedocs.io/en/latest/drizzlepac_api/astrodrizzle.html). Running `AstroDrizzle` using default parameter values is not recommended, as these defaults may not provide optimal science products. Users should also inspect the quality of the sky subtraction and cosmic ray rejection. For dithered data, users may experiment with the output `final_scale` and `final_pixfrac` parameters in the `final_drizzle` step. For more details, see the notebook in this repository to 'Optimize Image Sampling'." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# About this Notebook\n", + "## About this Notebook\n", + "\n", "\n", - " Author: C. Shanahan, STScI WFC3 Team \n", - " Updated: Jan. 20, 2022" + " Created: 14 Dec 2018; C. Shanahan & J. Mack\n", + " Updated: 31 Jul 2024; G. Anand & J. Mack \n", + "\n", + "\n", + "**Source:** GitHub [spacetelescope/hst_notebooks](https://github.com/spacetelescope/hst_notebooks)\n", + "\n", + "### Additional Resources \n", + "\n", + "Below are some additional resources that may be helpful. Please send any questions through the [HST Help Desk](https://stsci.service-now.com/hst), selecting the DrizzlePac category.\n", + "\n", + "- [ACS Website](https://www.stsci.edu/hst/instrumentation/acs)\n", + "- [ACS Instrument Handbook](https://hst-docs.stsci.edu/acsihb)\n", + "- [ACS Data Handbook](https://hst-docs.stsci.edu/acsdhb)\n", + "\n", + "- [WFC3 Website](https://www.stsci.edu/hst/instrumentation/wfc3)\n", + "- [WFC3 Instrument Handbook](https://hst-docs.stsci.edu/wfc3ihb)\n", + "- [WFC3 Data Handbook](https://hst-docs.stsci.edu/wfc3dhb)\n", + "\n", + "\n", + "### Citations \n", + "If you use Python packages such as `astropy`, `astroquery`, `drizzlepac`, `matplotlib`, or `numpy` for published research, please cite the authors. Follow these links for more information about citing various packages.\n", + "\n", + "* [Citing `astropy`](https://www.astropy.org/acknowledging.html)\n", + "* [Citing `astroquery`](https://github.com/astropy/astroquery/blob/main/astroquery/CITATION)\n", + "* [Citing `drizzlepac`](https://zenodo.org/records/3743274)\n", + "* [Citing `matplotlib`](https://matplotlib.org/stable/users/project/citing.html)\n", + "* [Citing `numpy`](https://numpy.org/citing-numpy/)\n", + "***\n", + "\n", + "[Top of Page](#top)\n", + "\"Space " ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -634,9 +1155,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.6" + "version": "3.11.7" } }, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 4 }