{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Aligning HST images to an absolute reference catalog\n", "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Note: The notebook in this repository 'Initializtion.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'.
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Note: This notebook is based on WFC3 ISR 2017-19: Aligning HST Images to Gaia: a Faster Mosaicking Workflow and contains a subset of the information/code found in the repository here. For more information, see the notebook in that repository titled 'Gaia_alignment.ipynb'.
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The alignment of HST exposures is a critical step in image stacking/combination performed by software such as `AstroDrizzle`. Generally, a relative alignment is performed that aligns one image (or multiple images) to another image which is designated as the reference image. This makes it so the images are aligned to each other, but the pointing error of the observatory can still cause the images to have incorrect absolute astrometry.\n", "\n", "When absolute astrometry is desired, the images can be aligned to an external catalog that is known to be on an absolute frame. In this example, we will provide a workflow to query catalogs such as SDSS and Gaia via the astroquery package, and then align the images to that catalog via TweakReg.\n", "\n", "For more information about TweakReg, see the other notebooks in this repository or the __[TweakReg Documentation](https://drizzlepac.readthedocs.io/en/deployment/tweakreg.html)__.\n", "\n", "For more information on Astroquery, see the other notebooks in this repository or the __[Astroquery Documentation](https://astroquery.readthedocs.io/en/latest/)__." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import astropy.units as u\n", "import glob\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import os\n", "\n", "from astropy.io import fits\n", "from astropy.table import Table\n", "from astropy.units import Quantity\n", "from astropy.coordinates import SkyCoord\n", "from astroquery.gaia import Gaia\n", "from astroquery.mast import Observations\n", "from astroquery.sdss import SDSS\n", "\n", "from ccdproc import ImageFileCollection\n", "from IPython.display import Image\n", "\n", "from drizzlepac import tweakreg\n", "from drizzlepac import astrodrizzle" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 1. Download the Data\n", "For this example, we will use WFC3/UVIS images of NGC 6791 from Visit 01 of proposal 12692." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Get the observation records\n", "obsTable = Observations.query_criteria(obs_id='ibwb01*', proposal_id=12692, obstype='all', filters='F606W')\n", "\n", "# Get the listing of data products\n", "products = Observations.get_product_list(obsTable)\n", "\n", "# Filter the products for exposures\n", "filtered_products = Observations.filter_products(products, productSubGroupDescription='FLC')\n", "\n", "# Show the table\n", "filtered_products" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Download all the images above\n", "Observations.download_products(filtered_products, mrp_only=False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# For convenience, move the products into the current directory.\n", "for flc in glob.glob('./mastDownload/HST/*/*flc.fits'):\n", " flc_name = os.path.split(flc)[-1]\n", " os.rename(flc, flc_name)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Inspect the image header\n", "\n", "The cell below shows how to query information from the image header using `ImageFileCollection` in `ccdproc`. \n", "We see that the 1st exposure is 30 seconds and the 2nd and 3rd exposures are 360 seconds. The 3rd exposure is dithered by ~82\" in the Y-direction which is approximately the width of one UVIS chip. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "collec = ImageFileCollection('./', glob_include=\"*flc.fits\", ext=0,\n", " keywords=[\"targname\", \"ra_targ\", \"dec_targ\", \"filter\", \"exptime\", \"postarg1\", \"postarg2\"])\n", "\n", "table = collec.summary\n", "table['exptime'].format = '7.1f'\n", "table['ra_targ'].format = '7.7f'\n", "table['dec_targ'].format = '7.7f'\n", "table['postarg1'].format = '7.2f'\n", "table['postarg2'].format = '7.2f'\n", "table" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 2. Querying catalogs\n", "\n", "Now that we have the images, we will download the reference catalogs from both SDSS and Gaia using `astroquery`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2a. Identify Coordinates\n", "We will first create a SkyCoord Object to point astroquery to where we are looking on the sky. Since our example uses data from NGC 6791, we will use the `ra_targ` and `dec_targ` keywords from the first image to get the coordinates of the object." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "RA = table['ra_targ'][0]\n", "Dec = table['dec_targ'][0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2b. SDSS Query\n", "We now give those values to an astropy `SkyCoord` object, which we will pass to the SDSS. Additionally, we use an astropy `Quantity` object to create a radius for the SDSS query. We set the radius to 6 arcminutes to comfortably cover the area of our images. For reference UVIS detector field of view is ~2.7'x2.7' and a y-dither of 82\" covers a total area on the sky of ~2.7'x4.1'." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "coord = SkyCoord(ra=RA, dec=Dec, unit=(u.deg, u.deg))\n", "radius = Quantity(6., u.arcmin)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we only need to perform the query via the `SDSS.query_region` method of `astroquery.sdss`. The `spectro=False` keyword argument means we want to exclude spectroscopic objects, as we are looking for objects to match with an image. \n", "\n", "In the fields parameter, we specify a list of fields we want returned by the query. In this case we only need the position, and maybe a magnitude 'g' if we want to cut very dim and/or bright objects out of the catalog, as those are likely measured poorly. Details on selecting objects by magnitude may be found in the original ['Gaia_alignment' notebook](https://github.com/spacetelescope/gaia_alignment). Many other fields are available in the SDSS query and are [documented here](http://cas.sdss.org/dr7/en/help/browser/description.asp?n=PhotoObj&t=V)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sdss_query = SDSS.query_region(coord, radius=radius, spectro=False, fields=['ra', 'dec', 'g'])\n", "sdss_query" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We then just need to save the catalog to use with TweakReg. Since the returned value of the query is an `astropy.Table`, saving the file is very straightforward:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sdss_query.write('sdss.cat', format='ascii.commented_header')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2c. Gaia Query\n", "Similarly to SDSS, we can query Gaia catalogs for our target via `astroquery.gaia`. We can use the same `coord` and `radius` from the SDSS query." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gaia_query = Gaia.query_object_async(coordinate=coord, radius=radius)\n", "gaia_query" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This query has returned very large number of columns. We want to pare down the catalog to make it easier to use with `TweakReg`. \n", "We can select only the useful columns via:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "reduced_query = gaia_query['ra', 'dec', 'phot_g_mean_mag']\n", "reduced_query" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we write this catalog to an ascii file for use with `TweakReg`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "reduced_query.write('gaia.cat', format='ascii.commented_header')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Running TweakReg\n", "With the catalogs downloaded and the headers populated, we simply need to run TweakReg with each catalog passed into the `refcat` parameter. The steps below compare the astrometric residuals obtained from aligning to each `refcat`. In each test case, we set `updatehdr` to False until we are satisfied with the alignment by inspecting both the shift file and the astrometric residual plots.\n", "\n", "### 3a. SDSS Alignment" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "refcat = 'sdss.cat'\n", "cw = 3.5 # Set to two times the FWHM of the PSF of the UVIS detector\n", "wcsname = 'SDSS' # Specify the WCS name for this alignment\n", "\n", "tweakreg.TweakReg('*flc.fits', # Pass input images\n", " updatehdr=False, # update header with new WCS solution\n", " imagefindcfg={'threshold': 500., 'conv_width': cw}, # Detection parameters, threshold varies for different data\n", " refcat=refcat, # Use user supplied catalog (Gaia)\n", " interactive=False,\n", " see2dplot=False,\n", " shiftfile=True, # Save out shift file (so we can look at shifts later)\n", " outshifts='SDSS_shifts.txt', # name of the shift file\n", " wcsname=wcsname, # Give our WCS a new name\n", " reusename=True,\n", " ylimit=1.5,\n", " fitgeometry='general') # Use the 6 parameter fit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can look at the shift file to see how well the fit did (or we could open the output png images for more information).\n", "\n", "The columns are:\n", "- Filename\n", "- X Shift [pixels]\n", "- Y Shift [pixels]\n", "- Rotation [degrees]\n", "- Scale\n", "- X RMS [pixels]\n", "- Y RMS [pixels]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for line in open('SDSS_shifts.txt').readlines():\n", " print(line)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Astrometric residual plots\n", "Image(filename='residuals_ibwb01xqq_flc.png',width=500, height=300)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Image(filename='residuals_ibwb01xrq_flc.png',width=500, height=300)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Image(filename='residuals_ibwb01xxq_flc.png',width=500, height=300)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see, the RMS is fairly large at about 0.5 pixels, which is not a great fit. This is likely because the SDSS astrometric precision is not high enough to get good HST alignment. One approach would be to align the first image to SDSS and then align the remaining HST images to one another. This would improve both the absolute and relative alignment of the individual frames.\n", "\n", "### 3b. Gaia Alignment" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "refcat = 'gaia.cat'\n", "cw = 3.5 # Set to two times the FWHM of the PSF.\n", "wcsname = 'Gaia' # Specify the WCS name for this alignment\n", "\n", "tweakreg.TweakReg('*flc.fits', # Pass input images\n", " updatehdr=False, # update header with new WCS solution\n", " imagefindcfg={'threshold':500.,'conv_width':cw}, # Detection parameters, threshold varies for different data\n", " refcat=refcat, # Use user supplied catalog (Gaia)\n", " interactive=False,\n", " see2dplot=False,\n", " shiftfile=True, # Save out shift file (so we can look at shifts later)\n", " outshifts='Gaia_shifts.txt', # name of the shift file\n", " wcsname=wcsname, # Give our WCS a new name\n", " reusename=True,\n", " sigma=2.3,\n", " ylimit=0.2,\n", " fitgeometry='general') # Use the 6 parameter fit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can similarly look at the shift file from alignment to the Gaia catalog:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for line in open('Gaia_shifts.txt').readlines():\n", " print(line)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Astrometric residual plots\n", "Image(filename='residuals_ibwb01xqq_flc.png',width=500, height=300)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Image(filename='residuals_ibwb01xrq_flc.png',width=500, height=300)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Image(filename='residuals_ibwb01xxq_flc.png',width=500, height=300)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected, the Gaia catalog does quite a bit better, with rms residuals less tha 0.05 pixels. \n", "\n", "To apply these transformations to the image, we simply need to run TweakReg the same as before, but set the parameter `updatehdr` equal to `True`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "refcat = 'gaia.cat'\n", "cw = 3.5 # Set to two times the FWHM of the PSF.\n", "wcsname = 'Gaia' # Specify the WCS name for this alignment\n", "\n", "tweakreg.TweakReg('*flc.fits', # Pass input images\n", " updatehdr=True, # update header with new WCS solution\n", " imagefindcfg={'threshold': 500., 'conv_width': cw}, # Detection parameters, threshold varies for different data\n", " refcat=refcat, # Use user supplied catalog (Gaia)\n", " interactive=False,\n", " see2dplot=False,\n", " shiftfile=True, # Save out shift file (so we can look at shifts later)\n", " outshifts='Gaia_shifts.txt', # name of the shift file\n", " wcsname=wcsname, # Give our WCS a new name\n", " reusename=True,\n", " sigma=2.3,\n", " fitgeometry='general') # Use the 6 parameter fit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Drizzle the Data\n", "\n", "While the three sets of FLC files are now aligned, we drizzle together only the two long exposures. \n", "\n", "When exposures are very different lengths, drizzling them together doesn't work well when 'EXP' weighting is used. For objects that saturate in the long exposures, the problem occurs at the boundary where you transition from only short exposure to short plus long. Here the pixels getting power from long exposure pixels are only getting power from pixels whose centers are outside the ring, and thus they are weighted lower than they would be if they were getting values from both inside and outside the ring. The result is a discontinuity in the PSF radial profile and a resulting flux which is too low in those boundary pixels. For photometry of saturated objects, the short exposures should be drizzled separately from the long exposures. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "astrodrizzle.AstroDrizzle('ibwb01x[rx]q_flc.fits', \n", " output='f606w',\n", " preserve=False,\n", " clean=True, \n", " build=False,\n", " context=False,\n", " skymethod='match',\n", " driz_sep_bits='64, 32',\n", " combine_type='minmed',\n", " final_bits='64, 32')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Display the combined science and weight images \n", "\n", "sci = fits.getdata('f606w_drc_sci.fits')\n", "wht = fits.getdata('f606w_drc_wht.fits')\n", "\n", "fig = plt.figure(figsize=(20, 20))\n", "ax1 = fig.add_subplot(1, 2, 1)\n", "ax2 = fig.add_subplot(1, 2, 2)\n", "\n", "ax1.imshow(sci, vmin=-0.05, vmax=0.4, cmap='Greys_r', origin='lower')\n", "ax2.imshow(wht, vmin=0, vmax=1000, cmap='Greys_r', origin='lower')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Conclusions\n", "\n", "Many other services have interfaces for querying catalogs which could also be used to align HST images. In general, Gaia works very well for HST due to it's high precision, but can have a low number of sources in some regions, especially at high galactic latitudes. Aligning images to an absolute frame provides an easy way to make data comparable across many epochs/detectors/observatories, and in many cases, makes the alignment much easier." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# About this Notebook\n", "\n", " Author: V. Bajaj, STScI WFC3 Team\n", " Updated: December 14, 2018" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.12" } }, "nbformat": 4, "nbformat_minor": 2 }