diff --git a/notebooks/DrizzlePac/sky_matching/MDRIZSKY_Values.png b/notebooks/DrizzlePac/sky_matching/MDRIZSKY_Values.png deleted file mode 100644 index 4b110fdb4..000000000 Binary files a/notebooks/DrizzlePac/sky_matching/MDRIZSKY_Values.png and /dev/null differ diff --git a/notebooks/DrizzlePac/sky_matching/drz.list b/notebooks/DrizzlePac/sky_matching/drz.list deleted file mode 100644 index bee0386d0..000000000 --- a/notebooks/DrizzlePac/sky_matching/drz.list +++ /dev/null @@ -1,9 +0,0 @@ -ibxl54030_drz.fits -ibxl51030_drz.fits -ibxl53030_drz.fits -ibxl55030_drz.fits -ibxl57030_drz.fits -ibxl50030_drz.fits -ibxl52030_drz.fits -ibxl56030_drz.fits -ibxl58030_drz.fits diff --git a/notebooks/DrizzlePac/sky_matching/requirements.txt b/notebooks/DrizzlePac/sky_matching/requirements.txt index 920a0035f..e7add88e2 100644 --- a/notebooks/DrizzlePac/sky_matching/requirements.txt +++ b/notebooks/DrizzlePac/sky_matching/requirements.txt @@ -1,9 +1,7 @@ -astropy==5.3.3 -astroquery==0.4.6 -ccdproc==2.4.0 -drizzlepac==3.5.1 -ipython==8.11.0 -matplotlib==3.7.0 -photutils==1.12.0 -pandas==1.5.3 -stwcs==1.7.2 +astropy>=5.3.3 +astroquery>=0.4.6 +drizzlepac>=3.5.1 +ipython>=8.11.0 +matplotlib>=3.7.0 +pandas>=1.5.3 +photutils==1.12.0 # The drizzlepac needs deprecated methods such as DAOGroup. \ No newline at end of file diff --git a/notebooks/DrizzlePac/sky_matching/sky_matching.ipynb b/notebooks/DrizzlePac/sky_matching/sky_matching.ipynb index e4703571a..5057ec5d2 100644 --- a/notebooks/DrizzlePac/sky_matching/sky_matching.ipynb +++ b/notebooks/DrizzlePac/sky_matching/sky_matching.ipynb @@ -4,14 +4,27 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Sky Matching " + "# Sky Matching for HST Mosaics \n", + "---" ] }, { "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 assumes you have created and activated a virtual environment using the requirements file in this notebook's repository. Please make sure you have read the contents of the README file before continuing the notebook. Note that the GIF file \"sky_matching_comparison.gif\" is one of the downloads needed for this notebook:
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Learning Goals:\n", + "In this tutorial we explore different options for handling the background sky with [`AstroDrizzle`](https://drizzlepac.readthedocs.io/en/latest/drizzlepac_api/astrodrizzle.html).
\n", + "By the end of this notebook you will:
\n", + "    • Download data with [`astroquery`](https://astroquery.readthedocs.io/en/latest/index.html)
\n", + "    • Align data with [`TweakReg`](https://drizzlepac.readthedocs.io/en/latest/user_reprocessing/tweakreg_api.html)
\n", + "    • Compare background sky options using the [`AstroDrizzle`](https://drizzlepac.readthedocs.io/en/latest/drizzlepac_api/astrodrizzle.html) parameter `skymethod`" ] }, { @@ -22,21 +35,29 @@ "\n", "[Introduction](#intro)
\n", "\n", - "[1. Retrieve Observations using Astroquery](#retrieve)
\n", - "- [Move files to the local directory](#move)
\n", - "- [Inspect the image headers](#inspect)
\n", - "\n", - "[2. Align the visit-level drizzled data](#align)
\n", - "\n", - "[3. Compare `skymethod` options in AstroDrizzle](#compare)
\n", - "- [`skymethod = 'localmin'`](#localmin)
\n", - "- [`skymethod = 'match'`](#match)
\n", - "- [`skymethod = 'globalmin+match'`](#globalminmatch)
\n", - "- [`skymethod = 'globalmin'`](#globalmin)
\n", - "\n", + "[1. Download the Observations with `astroquery`](#download)
\n", + "    [1.1 Check image header data](#check_keywords)
\n", + "    [1.2 Inspect the alignment](#check_wcs)
\n", + "[2. Align the visit-level drizzled data with `TweakReg`](#align)
\n", + "    [2.1 Create a catalog of Gaia DR3 sources](#cat)
\n", + "    [2.2 Create a catalog of Gaia DR3 sources with Proper Motion Data](#cat_pm)
\n", + "    [2.3 Run `Tweakreg`](#tweak)
\n", + "    [2.4 Inspect the shift file and fit quality](#fit_quality)
\n", + "    [2.5 Overplot matched sources and inspect fit residuals](#overplot)
\n", + "    [2.6 Rerun `TweakReg` and update the header WCS](#updatehdr)
\n", + "    [2.7 Run `TweakBack` to propogate the WCS to the FLT files](#tweakback)
\n", + "[3. Compare `skymethod` options in `AstroDrizzle`](#compare)
\n", + "    [3.1 `skymethod = 'localmin'`](#localmin)
\n", + "    [3.2 `skymethod = 'match'`](#match)
\n", + "    [3.3 `skymethod = 'globalmin+match'`](#globalminmatch)
\n", + "    [3.4 `skymethod = 'globalmin'`](#globalmin)
\n", "[4. Compare the MDRIZSKY values for each method](#mdrizsky)
\n", + "[5. Display the 'sky matched' science mosaic and weight image](#display)
\n", + "[6. Conclusion](#conclusion)
\n", "\n", - "[5. Display the full 'matched' mosaic](#display)
" + "[Additional Resources](#resources)
\n", + "[About this notebook](#about)
\n", + "[Citations](#citations)" ] }, { @@ -45,17 +66,17 @@ "source": [ "## Introduction \n", "\n", - "When a user decided to mosaic several images it can be difficult to decide how to relate the background of one image to the others. AstroDrizzle has the ability to compute and then either subtract or equalize the background in input images in order to match the sky between those images. The function allows users to select the algorithm for the sky subtraction using the `skymethod` parameter. \n", + "When creating an image mosaic, `AstroDrizzle` has the ability to compute the sky and then either subtract or equalize the background in input images. Users may select the algorithm used for the sky subtraction via the `skymethod` parameter. \n", "\n", - "There are four methods available in sky matching: 'localmin', 'match', 'globalmin', 'globalmin+match'.\n", + "There are four methods available in sky matching: `localmin`, `match`, `globalmin`, and `globalmin+match`.\n", "\n", - "By applying `drizzlepac.sky.sky()`, or by applying the `skymethod` parameter in the call to `drizzlepac.astrodrizzle.AstroDrizzle()`, AstroDrizzle will update the keyword `MDRIZSKY` in the headers of the input files. \n", + "By applying `drizzlepac.sky.sky()`, or using the `skymethod` parameter in the call to `drizzlepac.astrodrizzle.AstroDrizzle()`, AstroDrizzle will update the keyword `MDRIZSKY` in the headers of the input files but it will not change the science data. \n", "\n", - "The **recommended** use for images with sparse fields that have few astronomical sources is `skymethod = localmin` and for images with complicated backgrounds, such as nebulae and large host galaxies, `skymethod = match` is recommended.\n", + "For images of sparse fields with few astronomical sources, the default `skymethod = 'localmin'` may be used, although this method can slightly oversubtract the background. For images with complicated backgrounds, such as nebulae and large host galaxies, `skymethod = 'match'` is recommended.\n", "\n", - "For more information on the specifics of this function please refer to the documentation [here](https://drizzlepac.readthedocs.io/en/deployment/sky.html#drizzlepac.sky.sky).\n", + "For more information on the specifics of this function, please refer to the documentation [here](https://drizzlepac.readthedocs.io/en/latest/drizzlepac_api/sky.html#drizzlepac.sky.sky)\n", "\n", - "Below, each of the four methods will be demonstrated using a single example dataset, and differences between the methods will be highlighted. " + "Below, each of the four methods is demonstrated using a single example dataset, and differences between the methods is highlighted. " ] }, { @@ -65,66 +86,213 @@ "outputs": [], "source": [ "# All imports needed through out this notebook are included at the beginning. \n", - "import glob\n", - "from shutil import copyfile\n", + "%matplotlib inline\n", + "%config InlineBackend.figure_format = 'retina'\n", "\n", - "from astropy.io import fits\n", - "from astroquery.mast import Observations\n", - "from ccdproc import ImageFileCollection\n", + "from collections import defaultdict\n", + "from IPython.display import clear_output \n", + "import glob\n", + "import os\n", + "import shutil \n", "import matplotlib.pyplot as plt\n", - "import pandas as pd\n", + "import matplotlib.image as mpimg\n", + "import numpy as np\n", + "import pandas\n", + "\n", + "from astropy.coordinates import SkyCoord\n", + "from astropy.io import ascii, fits\n", + "from astropy.table import Table\n", + "from astropy.units import Quantity\n", + "import astropy.units as u\n", + "from astroquery.gaia import Gaia\n", + "from astroquery.mast import Observations\n", + "from drizzlepac import astrodrizzle, tweakback, tweakreg \n", + "from drizzlepac.haputils.astrometric_utils import create_astrometric_catalog\n", + "\n", "\n", - "import drizzlepac\n", - "from drizzlepac import tweakreg\n", - "from drizzlepac import tweakback" + "Gaia.MAIN_GAIA_TABLE = 'gaiadr3.gaia_source' # Change if different data release is desired\n", + "Gaia.ROW_LIMIT = 100000" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## 1. Retrieve Observations using Astroquery " + "## 1. Download the Observations with `astroquery` \n", + "\n", + "MAST queries may be done using `query_criteria`, where we specify:
\n", + "\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", + "WFC3/IR observations of the Horsehead Nebula in the F160W filter obtained in HST proposal\n", + "program [12812](http://www.stsci.edu/cgi-bin/get-proposal-info?id=12812&observatory=HST)\n", + "will be used for this demonstration. \n", + "\n", + "Nine visits were acquired in a 3x3 mosaic pattern on the sky, with two dither positions per visit in two IR filters. [High level science products](https://archive.stsci.edu/prepds/heritage/horsehead/readme_HLSP_v3.txt) for these datasets were delivered to MAST in 2013, and this notebook is based on that user tutorial but has been updated to align these data to Gaia.\n", + "\n", + "The 18 FLT images ibxl5*_flt.fits have been processed by the HST WFC3 pipeline (calwf3), which includes bias subtraction, dark current correction, cosmic-ray rejection, and flatfielding. The 9 DRZ files ibxl5*_drz.fits have been processed with `AstroDrizzle` to remove distortion and to combine the 2 dithered FLT frames by filter for each vist.\n", + "\n", + "
THIS IS A LARGE DOWNLOAD (~400 MB). Depending on your connection speed, the next cell may take a few minutes to execute.
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "obs_ids = ['ibxl5*']\n", + "props = ['12812']\n", + "filts = ['F160W']\n", + "\n", + "obsTable = Observations.query_criteria(obs_id=obs_ids, proposal_id=props, filters=filts)\n", + "products = Observations.get_product_list(obsTable)\n", + "\n", + "data_prod = ['FLT', 'DRZ'] # ['FLC','FLT','DRC','DRZ'] \n", + "data_type = ['CALWF3'] # ['CALACS','CALWF3','CALWP2','HAP-SVM'] \n", + "\n", + "Observations.download_products(products, download_dir='./science',\n", + " productSubGroupDescription=data_prod, \n", + " project=data_type)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Move the files to the local working directory" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "files = glob.glob(os.path.join(os.curdir, 'science', 'mastDownload', 'HST', '*', '*fits'))\n", + "for im in files:\n", + " root = os.path.basename(im)\n", + " os.rename(im, './' + root)\n", + " \n", + "if os.path.exists('./science'):\n", + " shutil.rmtree('science/')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "WFC3/IR observations of the Horsehead Nebula obtained in HST proposal 12812 (PI: Levay) will be used for this demonstration. The images have been processed by the HST WFC3 pipeline (calwf3), which includes bias subtraction, dark current correction, cosmic-ray rejection, and flatfielding.\n", + "### 1.1 Check image header data \n", "\n", - "The alignment specifics for this example are based on the example found [here](https://archive.stsci.edu/prepds/heritage/horsehead/readme_HLSP_v3.txt). We will provide explanation about the process as we go. \n", + "Here we will look at important keywords in the image headers." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "files = sorted(glob.glob('*fl?.fits'))\n", + "keywords_ext0 = [\"ROOTNAME\", \"ASN_ID\", \"TARGNAME\", \"DETECTOR\", \"FILTER\", \"EXPTIME\", \n", + " \"RA_TARG\", \"DEC_TARG\", \"POSTARG1\", \"POSTARG2\", \"DATE-OBS\"]\n", + "keywords_ext1 = [\"ORIENTAT\"]\n", + "data = []\n", + "\n", + "for file in files:\n", + " path_data = []\n", + " for keyword in keywords_ext0:\n", + " path_data.append(fits.getval(file, keyword, ext=0))\n", + " for keyword in keywords_ext1:\n", + " path_data.append(fits.getval(file, keyword, ext=1))\n", + " data.append(path_data)\n", + " \n", + "keywords = keywords_ext0 + keywords_ext1\n", + "table = Table(np.array(data), names=keywords, dtype=['str', 'str', 'str', 'str', 'str', 'f8', 'f8', 'f8', 'f8', 'f8', 'str', 'f8'])\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.show_in_notebook()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1.2 Inspect the Alignment \n", "\n", - "The code used to retrieve the data from astroquery is below." + "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": { - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ - "# Retrieve the observation information.\n", - "obs_table = Observations.query_criteria(proposal_id='12812', filters=[\"F160W\"], project='HST')\n", - "obs_table.show_in_notebook()\n", + "ext_0_kws = ['DETECTOR']\n", + "ext_1_kws = ['WCSNAME', 'NMATCHES', 'RMS_RA', 'RMS_DEC']\n", "\n", - "# Figure out the downloadable products:\n", - "dataProductsByObservation = Observations.get_product_list(obs_table)\n", - "dataProductsByObservation.show_in_notebook()\n", + "det_scale = {'IR': 0.1283, 'UVIS': 0.0396, 'WFC': 0.05} # plate scale (arcsec/pixel)\n", "\n", - "# Download the flc.fits files: \n", - "obsids = obs_table['obsid']\n", + "format_dict = {}\n", + "col_dict = defaultdict(list)\n", "\n", - "# Download the data - uncomment to download: \n", - "Observations.download_products(obsids, mrp_only=False, productSubGroupDescription=['FLT', 'DRZ'])" + "for f in sorted(glob.glob('*dr?.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(hdr1[kw]/1000./det_scale[hdr0['DETECTOR']], 2) # convert RMS from mas to pixels\n", + " col_dict[f'{kw}_pix'].append(val)\n", + "\n", + "wcstable = Table(col_dict)\n", + "wcstable.show_in_notebook()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Move files to the local directory " + "
Note that there are different WCS solutions for each visit, with Gaia eDR3 as the reference catalog for all but two which were fit to GSC v2.4.2 and which have a much larger fit rms values (>0.5 pixels). Since the WCS solutions are inconsistent for this target, we wish to realign the data to use a common reference catalog.
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Align the visit-level drizzled data with `TweakReg` \n", + "\n", + "Here we will use `TweakReg` to align the DRZ files to Gaia DR3 and then use `TweakBack` to propagate those solutions back to the FLT image headers prior to combining with `AstroDrizzle`. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2.1 Create a catalog of Gaia DR3 sources \n", + "\n", + "This method uses the RA/Dec of the first image and a radius of 5'." ] }, { @@ -133,28 +301,61 @@ "metadata": {}, "outputs": [], "source": [ - "flt_files = glob.glob('mastDownload/*/*/*flt.fits')\n", + "RA = table['RA_TARG'][0]\n", + "Dec = table['DEC_TARG'][0]\n", "\n", - "for file in flt_files:\n", - " with fits.open(file) as im:\n", - " rootname = im[0].header['ROOTNAME']\n", - " dst = rootname + \"_flt.fits\"\n", - " copyfile(file, dst)\n", + "coord = SkyCoord(ra=RA, dec=Dec, unit=(u.deg, u.deg))\n", + "radius = Quantity(5., u.arcmin)\n", "\n", - "drz_files = glob.glob('mastDownload/*/*/*drz.fits')\n", + "gaia_query = Gaia.query_object_async(coordinate=coord, radius=radius)\n", + "gaia_query" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "reduced_query = gaia_query['ra', 'dec', 'phot_g_mean_mag']\n", + "reduced_query.write('gaia_no_pm.cat', format='ascii.commented_header', overwrite=True)\n", + "reduced_query" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2.2 Create a catalog of Gaia DR3 sources with Proper Motion Data \n", "\n", - "for file in drz_files:\n", - " with fits.open(file) as im:\n", - " rootname = im[0].header['ROOTNAME']\n", - " dst = rootname + \"_drz.fits\"\n", - " copyfile(file, dst)" + "This method uses the image FLT footprints and gives 161 sources, compared to 183 with the prior method." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pm_cat = create_astrometric_catalog(sorted(glob.glob('*flt.fits')))\n", + "pm_cat.write('gaia_pm.cat', overwrite=True, format='ascii.no_header')\n", + "len(pm_cat)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Inspect the image headers " + "### 2.3 Run `TweakReg` " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we run `TweakReg` on the visit-level drizzled (DRZ) images and align to the Gaia catalog with proper motion data. \n", + "\n", + "Because the fit RMS values for the MAST products were large for some visits, we allow for a larger than usual search radius of 1\". We also set the `conv_width` value slightly higher than the recommended value of 2.5 for WFC3/IR data in order to use barely resolved sources for alignment." ] }, { @@ -163,15 +364,42 @@ "metadata": {}, "outputs": [], "source": [ - "# FLT images\n", - "collect_flt = ImageFileCollection('./', glob_include=\"*flt.fits\", ext=0,\n", - " keywords=[\"asn_id\", \"detector\", \"filter\", \"exptime\", \"postarg1\", \"postarg2\"])\n", - " \n", - "flt_table = collect_flt.summary\n", - "flt_table['exptime'].format = '7.1f'\n", - "flt_table['postarg1'].format = '7.2f'\n", - "flt_table['postarg2'].format = '7.2f'\n", - "flt_table" + "refcat = 'gaia_pm.cat' # Use the catalog with proper motion data\n", + "\n", + "drz_files = sorted(glob.glob('*drz.fits'))\n", + "\n", + "tweakreg.TweakReg(drz_files, \n", + " imagefindcfg={'threshold': 4, 'conv_width': 4.5}, \n", + " minobj=3,\n", + " shiftfile=True, \n", + " outshifts='shift_drz.txt',\n", + " refcat=refcat,\n", + " searchrad=1,\n", + " ylimit=0.4, \n", + " nclip=1,\n", + " updatehdr=False, # change later when you verify the alignment works\n", + " interactive=False)\n", + "# clear_output()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# If the alignment is unsuccessful, stop the notebook\n", + "with open('shift_drz.txt', 'r') as shift:\n", + " for line_number, line in enumerate(shift, start=1):\n", + " if \"nan\" in line:\n", + " raise ValueError(f'nan found in line {line_number} in shift file')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2.4 Inspect the shift file and fit quality " ] }, { @@ -180,61 +408,153 @@ "metadata": {}, "outputs": [], "source": [ - "# DRZ images\n", - "collect_drz = ImageFileCollection('./', glob_include=\"*drz.fits\", ext=0,\n", - " keywords=[\"asn_id\", \"detector\", \"filter\", \"exptime\"])\n", + "# Read the shift file just created by tweakreg\n", + "# There are 7 columns including: filename, x-shift, y-shift, rotation, scale, x-RMS, and y-RMS\n", + "shift_table = Table.read('shift_drz.txt',\n", + " format='ascii.no_header', \n", + " names=['file', 'dx', 'dy', 'rot', 'scale', 'xrms', 'yrms'])\n", + "\n", + "# Define the format for each column (excluding 'file').\n", + "formats = ['.2f', '.2f', '.3f', '.5f', '.2f', '.2f']\n", + "\n", + "# Iterate over the columns 'dx', 'dy', 'rot', 'scale', 'xrms', 'yrms'\n", + "for i, col in enumerate(shift_table.colnames[1:]):\n", + " # Apply the format to the current column \n", + " shift_table[col].format = formats[i]\n", " \n", - "drz_table = collect_drz.summary\n", - "drz_table['exptime'].format = '7.1f'\n", - "drz_table" + "# Display the table in the notebook\n", + "shift_table" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## 2. Align the visit-level drizzled data " + "Note that there are large residual rotation and scale terms in the shift file for several visits in the MAST data products. These will be corrected when we run `TweakReg` an additional time and update the WCS in Section 2.6 below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "match_files = sorted(glob.glob('*_drz_catalog_fit.match'))\n", + "for f in match_files:\n", + " input = ascii.read(f)\n", + " print(f'Number of matches for {f} {len(input)}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Exposures obtained within a single HST visit tend to be aligned very well since they use the same guide stars. Thus, we will use TweakReg to align the DRZ files and then use TweakBack to propagate those solutions back to the FLT image headers prior to combining with AstroDrizzle. Making use of the parameter `expand_refcat`, TweakReg will build up an expanded reference catalog on the sky to be used for alignment. For this dataset, we obtain the best results when we tell TweakReg to align the DRZ files in a specific order, and this is achieved by giving a list of files as input. More details on alignment of HST mosaics can be found in the notebook 'mosaics.ipynb' in this repository.\n", - " \n", - " drz.list\n", - " ibxl54030_drz.fits\t<-- Note that tile 54 is the reference and is listed first\n", - " ibxl51030_drz.fits \n", - " ibxl53030_drz.fits\n", - " ibxl55030_drz.fits\n", - " ibxl57030_drz.fits\n", - " ibxl50030_drz.fits\n", - " ibxl52030_drz.fits\n", - " ibxl56030_drz.fits\n", - " ibxl58030_drz.fits\n", - " ____ ____ ____\n", - " | | | | <-- The 4 tiles with the most overlap are aligned first\n", - " | | 55 | |\n", - " |____|____|____|\n", - " | | | |\n", - " | 57 | 54 | 51 |\n", - " |____|____|____|\n", - " | | | |\n", - " | | 53 | |\n", - " |____|____|____|\n", - " \n", - " ____ ____ ____\n", - " | | | | <-- Then the corner tiles are added \n", - " | 58 | | 52 |\n", - " |____|____|____|\n", - " | | | |\n", - " | | | |\n", - " |____|____|____|\n", - " | | | |\n", - " | 56 | | 50 |\n", - " |____|____|____|\n", - "\n", - "This allows TweakReg to start with the tiles with the most overlap first and then build upon that, expanding the reference catalog as each new tile is aligned. " + "Next we compare with the MAST WCS solutions, number of matches and fit RMS values. Since there are only 5 Gaia sources in Visits 53 and 56, these did not have a successful Gaia fit during MAST processing and instead were aligned to the next catalog in the priority list, GSC v2.4.2. (Currently the minimum number of matches for a successful fit is 6, but this will updated to 10 in summer 2024. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "wcstable" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2.5 Overplot matched sources and inspect fit residuals " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we overplot the HST sources which were successfully matched with Gaia eDR3 and we look at the astrometric fit residual PNG plots.\n", + "\n", + "While we inspect only two visits (52 and 53) at a time, and additional visits may be uncommented in the cell below. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define the rootnames for the two FITS files to be compared\n", + "# Uncomment the pairs you want to use\n", + "\n", + "# rootname_A = 'ibxl50030'\n", + "# rootname_B = 'ibxl51030'\n", + "rootname_A = 'ibxl52030'\n", + "rootname_B = 'ibxl53030'\n", + "# rootname_A = 'ibxl54030'\n", + "# rootname_B = 'ibxl55030'\n", + "# rootname_A = 'ibxl56030'\n", + "# rootname_B = 'ibxl57030'\n", + "# rootname_A = 'ibxl57030'\n", + "# rootname_B = 'ibxl58030'\n", + "\n", + "# Create subplots with 1 row and 2 columns\n", + "fig, [ax1, ax2] = plt.subplots(1, 2, figsize=(13, 10))\n", + "\n", + "# Open the FITS files and extract the data from the 'SCI' extension\n", + "data_a = fits.open(rootname_A+'_drz.fits')['SCI', 1].data\n", + "data_b = fits.open(rootname_B+'_drz.fits')['SCI', 1].data\n", + "\n", + "# Display the images in grayscale with a stretch from 0 to 2 \n", + "ax1.imshow(data_a, cmap='Greys', origin='lower', vmin=0, vmax=2)\n", + "ax2.imshow(data_b, cmap='Greys', origin='lower', vmin=0, vmax=2)\n", + "\n", + "# Read the matching catalog files for both rootnames\n", + "match_tab_a = ascii.read(rootname_A+'_drz_catalog_fit.match')\n", + "match_tab_b = ascii.read(rootname_B+'_drz_catalog_fit.match')\n", + "\n", + "# Extract x and y coordinates from the matching catalogs\n", + "x_coord_a, y_coord_a = match_tab_a['col11'], match_tab_a['col12']\n", + "x_coord_b, y_coord_b = match_tab_b['col11'], match_tab_b['col12']\n", + "\n", + "# Plot the coordinates on the images with red circles\n", + "ax1.scatter(x_coord_a, y_coord_a, s=30, edgecolor='r', facecolor='None')\n", + "ax2.scatter(x_coord_b, y_coord_b, s=30, edgecolor='r', facecolor='None')\n", + "\n", + "# Set the titles for the subplots, including the number of matches\n", + "ax1.set_title(rootname_A+f' N = {len(match_tab_a)} Gaia Matches', fontsize=20)\n", + "ax2.set_title(rootname_B+f' N = {len(match_tab_b)} Gaia Matches', fontsize=20)\n", + "fig.tight_layout()\n", + "\n", + "# Load, display, and inspect fit residual PNG files\n", + "img_A = mpimg.imread(f'residuals_{rootname_A}_drz.png')\n", + "img_B = mpimg.imread(f'residuals_{rootname_B}_drz.png')\n", + "\n", + "# Create subplots for the residual 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", + "\n", + "# Remove unnecessary axis from the residual images\n", + "axs[0].axis('off') \n", + "axs[1].axis('off')\n", + "\n", + "# Adjust layout to minimize margins and display plots\n", + "fig.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2.6 Rerun `TweakReg` and update the header WCS " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Once we are happy with the alignment, we run `TweakReg` again with the same parameters but with the `updatehdr=True`." ] }, { @@ -245,26 +565,28 @@ }, "outputs": [], "source": [ - "tweakreg.TweakReg('@drz.list', \n", - " imagefindcfg={'threshold': 5, 'conv_width': 4.5, 'peakmin': 1}, \n", - " minobj=5, \n", - " shiftfile=True, \n", - " expand_refcat=True, \n", - " enforce_user_order=True, \n", - " outshifts='shift4_drc.txt',\n", + "refcat = 'gaia_pm.cat' # Use the catalog with proper motion data\n", + "\n", + "tweakreg.TweakReg(drz_files, \n", + " imagefindcfg={'threshold': 4, 'conv_width': 4.5}, \n", + " minobj=4,\n", + " shiftfile=False, \n", + " refcat=refcat,\n", " searchrad=1,\n", - " ylimit=0.5, \n", - " tolerance=1,\n", + " ylimit=0.4, \n", " nclip=1,\n", - " updatehdr=True,\n", - " interactive=False)" + " updatehdr=True, # update header\n", + " interactive=False)\n", + "clear_output()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Tweakback is then run on the aligned DRZ files to propogate the updated WCS information back to the FLT files." + "### 2.7 Run `TweakBack` to propogate the WCS to the FLT files \n", + "\n", + "Finally, we run `Tweakback` on the aligned DRZ files to propogate the updated WCS information back to the FLT files." ] }, { @@ -275,37 +597,46 @@ }, "outputs": [], "source": [ - "drz_files = glob.glob('i*0_drz.fits')\n", - "\n", "for f in drz_files:\n", " tb_input = f+'[sci,1]'\n", - " hdu = fits.open(f)\n", - " tweakback.apply_tweak(tb_input, orig_wcs_name=hdu[1].header['WCSNAME'])" + " with fits.open(f) as hdu:\n", + " tweakback.apply_tweak(tb_input, orig_wcs_name=hdu[1].header['WCSNAME'])\n", + " \n", + "clear_output()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## 3. Compare `skymethod` options in AstroDrizzle " + "## 3. Compare `skymethod` options in AstroDrizzle \n", + "\n", + "Now that the FLT files contain the updated WCS solutions, we explore different algorithms for estimating the sky." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### `skymethod = 'localmin'` " + "### 3.1 `skymethod = 'localmin'` " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "When using AstroDrizzle to compute the sky in each frame, 'localmin' will compute a common sky for all members of an exposure, as described [here](https://drizzlepac.readthedocs.io/en/deployment/sky.html#drizzlepac.sky.sky). This algorithm is recommended when images are dominated by blank sky instead of extended, diffuse sources.\n", + "When using `AstroDrizzle` to compute the sky in each frame, 'localmin' will compute a common sky value for all chips of a given exposure, using the minimum sky value from all chips. This process is repeated for each input image. This algorithm is recommended when images are dominated by blank sky instead of extended, diffuse sources.\n", "\n", - "For a typical use, it will compute sky values for each chip/image extension (marked for sky subtraction in the input parameter) in an input image, and it will subtract the previously found minimum sky value from all chips in that image. This process is repeated for each input image.\n", + "See [readthedocs](https://drizzlepac.readthedocs.io/en/deployment/sky.html#drizzlepac.sky.sky) for more details on sky subtraction options." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the command below, the aligned FLT frames are sky subtracted and drizzled together. \n", "\n", - "In the command below, the aligned FLT frames are sky subtracted and drizzled together. Because the WFC3/IR data products are already cleaned of cosmic rays during calwf3 processing, cosmic-ray rejection is turned off in AstroDrizzle by setting the parameters `driz_separate`, `median`, `blot`, and `driz_cr` to 'False'." + "Because the WFC3/IR data products are already cleaned of cosmic rays during calwf3 processing, cosmic-ray rejection is turned off in AstroDrizzle by setting the parameters `driz_separate`, `median`, `blot`, and `driz_cr` to 'False'. Note that `final_bits=16` means only the stable hot pixels are treated as good. " ] }, { @@ -316,24 +647,24 @@ }, "outputs": [], "source": [ - "drizzlepac.astrodrizzle.AstroDrizzle('*flt.fits', \n", - " output='f160w_localmin',\n", - " preserve=False, \n", - " skymethod='localmin', \n", - " driz_separate=False, \n", - " median=False, \n", - " blot=False, \n", - " driz_cr=False, \n", - " final_bits='64',\n", - " final_wcs=True, \n", - " final_rot=257.)" + "sky = 'localmin'\n", + "astrodrizzle.AstroDrizzle('*flt.fits', \n", + " output='f160w_localmin',\n", + " preserve=False, \n", + " context=False,\n", + " skymethod=sky, \n", + " driz_separate=False, median=False, blot=False, driz_cr=False, # CR-rej turned off\n", + " final_bits='16',\n", + " final_wcs=True, \n", + " final_rot=257.)\n", + "clear_output()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### `skymethod = 'match'` " + "### 3.2 `skymethod = 'match'` " ] }, { @@ -351,24 +682,24 @@ "metadata": {}, "outputs": [], "source": [ - "drizzlepac.astrodrizzle.AstroDrizzle('*flt.fits', \n", - " output='f160w_match',\n", - " preserve=False, \n", - " skymethod='match',\n", - " driz_separate=False, \n", - " median=False, \n", - " blot=False, \n", - " driz_cr=False, \n", - " final_bits='64',\n", - " final_wcs=True, \n", - " final_rot=257.)" + "sky = 'match'\n", + "astrodrizzle.AstroDrizzle('*flt.fits', \n", + " output='f160w_match',\n", + " preserve=False, \n", + " context=False,\n", + " skymethod=sky, \n", + " driz_separate=False, median=False, blot=False, driz_cr=False, # CRREJ=None\n", + " final_bits='16',\n", + " final_wcs=True, \n", + " final_rot=257.)\n", + "clear_output()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### `skymethod = 'globalmin+match'` " + "### 3.3 `skymethod = 'globalmin+match'` " ] }, { @@ -384,24 +715,24 @@ "metadata": {}, "outputs": [], "source": [ - "drizzlepac.astrodrizzle.AstroDrizzle('*flt.fits', \n", - " output='f160w_globalmin_match',\n", - " preserve=False, \n", - " skymethod='globalmin+match',\n", - " driz_separate=False, \n", - " median=False, \n", - " blot=False, \n", - " driz_cr=False, \n", - " final_bits='64',\n", - " final_wcs=True, \n", - " final_rot=257.)" + "sky = 'globalmin+match'\n", + "astrodrizzle.AstroDrizzle('*flt.fits', \n", + " output='f160w_globalmin_match',\n", + " preserve=False, \n", + " context=False,\n", + " skymethod=sky, \n", + " driz_separate=False, median=False, blot=False, driz_cr=False, # CRREJ=None\n", + " final_bits='16',\n", + " final_wcs=True, \n", + " final_rot=257.)\n", + "clear_output()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### `skymethod = 'globalmin'` " + "### 3.4 `skymethod = 'globalmin'` " ] }, { @@ -416,20 +747,22 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ - "drizzlepac.astrodrizzle.AstroDrizzle('*flt.fits', \n", - " output='f160w_globalmin', \n", - " preserve=False,\n", - " skymethod='globalmin', \n", - " driz_separate=False, \n", - " median=False, \n", - " blot=False, \n", - " driz_cr=False, \n", - " final_bits='64',\n", - " final_wcs=True, \n", - " final_rot=257.)" + "sky = 'globalmin'\n", + "astrodrizzle.AstroDrizzle('*flt.fits', \n", + " output='f160w_globalmin',\n", + " preserve=False, \n", + " context=False,\n", + " skymethod=sky, \n", + " driz_separate=False, median=False, blot=False, driz_cr=False, # CRREJ=None\n", + " final_bits='16',\n", + " final_wcs=True, \n", + " final_rot=257.)\n", + "clear_output()" ] }, { @@ -443,16 +776,16 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Below we provide a gif comparing the upper portion of the final drizzled image. We cycle through three of the versions that use different `skymethod` algorithms: \n", + "Below we provide a gif comparing the upper portion of the final drizzled image. We cycle through three drizzled images produced using different `skymethod` algorithms: \n", "\n", - "![SegmentLocal](labeled_local_globalmatch_match.gif)" + "![The top of the horsehead nebula is shown in grayscale with a stretch of about -0.4 to 2 and the animated gif transitions between displaying the \"localmin\", \"match\", and \"globalmin+match\" sky methods from astrodrizzle. Localmin looks the worst with unmatched backgrounds in the top left and right corners where the outline of separate pointings (dithers) are apparent.](labeled_local_globalmatch_match.gif)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Below is the code used to compare results from the various skymatching methods. It's best practice to use more compact for loops when possible in situations like this, but we wanted to write out each file to ensure we are gathering the proper `MDRIZSKY` values for each `skymethod` run. We also display a table (Pandas DataFrame) of the different sky matching `MDRIZSKY` values for comparison." + "Next we print the sky values computed for each image using the four different methods." ] }, { @@ -461,24 +794,11 @@ "metadata": {}, "outputs": [], "source": [ - "with fits.open(\"f160w_globalmin_drz_sci.fits\") as im:\n", - " globalmin_data = im[1].data\n", - "\n", - "with fits.open(\"f160w_globalmin_match_drz_sci.fits\") as im:\n", - " globalmin_match_data = im[1].data\n", - "\n", - "with fits.open(\"f160w_match_drz_sci.fits\") as im:\n", - " match_data = im[1].data\n", - "\n", - "with fits.open(\"f160w_localmin_drz_sci.fits\") as im:\n", - " localmin_data = im[1].data\n", - "\n", - "mdrizsky_val = pd.DataFrame(\n", - " {'Rootname': globalmin_data[\"ROOTNAME\"],\n", - " 'globalmin': globalmin_data[\"MDRIZSKY\"],\n", - " 'globalmin_match': globalmin_match_data[\"MDRIZSKY\"],\n", - " 'match': match_data[\"MDRIZSKY\"],\n", - " 'local': localmin_data['MDRIZSKY']})\n", + "mdrizsky_val = pandas.DataFrame({'rootname': fits.getdata('f160w_globalmin_drz_sci.fits', 1)['rootname'],\n", + " 'local': fits.getdata('f160w_localmin_drz_sci.fits', 1)['mdrizsky'],\n", + " 'globalmin': fits.getdata('f160w_globalmin_drz_sci.fits', 1)['mdrizsky'],\n", + " 'globalmin_match': fits.getdata('f160w_globalmin_match_drz_sci.fits', 1)['mdrizsky'],\n", + " 'match': fits.getdata('f160w_match_drz_sci.fits', 1)['mdrizsky']})\n", "mdrizsky_val" ] }, @@ -486,7 +806,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "These computed sky values can be visualized in the plot below. To reiterate, the MDRIZSKY keyword reports the value subtracted from each image, and not the sky level itself. Thus the values for `skymethod='match'` are close to zero. We also note that varying background levels across the individual tiles result in inaccurate sky background determination when `skymethod='match'` and thus a mismatched sky in the final mosaic." + "These computed sky values can be visualized in the plot below. To reiterate, the MDRIZSKY keyword reports the value subtracted from each FLC image prior to drizzling, and not the sky level itself. Thus the values for `skymethod='match'` are close to zero. \n", + "\n", + "We also note that varying background levels across the individual tiles result in inaccurate sky background determination when `skymethod='localmin'` and thus a mismatched sky in the final mosaic." ] }, { @@ -502,14 +824,14 @@ "local = list(mdrizsky_val['local'])\n", "\n", "# Plotting code: \n", - "fig = plt.figure(figsize=[10, 10])\n", + "fig = plt.figure(figsize=[7, 7])\n", "plt.scatter(index, globalmin_match, color='magenta', label='Globalmin + Match')\n", "plt.scatter(index, match, color='navy', label='Match')\n", "plt.scatter(index, local, color='olive', label='Local')\n", "plt.scatter(index, globalmin, color='orange', label='Globalmin')\n", "plt.xlabel('Individual Images')\n", "plt.ylabel('MDRIZSKY Value')\n", - "plt.legend(loc=\"center left\")\n", + "plt.legend(loc='best')\n", "plt.xticks(index)\n", "plt.tight_layout()\n", "plt.show()" @@ -519,7 +841,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## 5. Display the full 'matched' mosaic " + "## 5. Display the 'sky matched' science mosaic and weight image \n", + "\n", + "Finally, we display the science and weight images for the combined mosaic." ] }, { @@ -529,30 +853,89 @@ "outputs": [], "source": [ "sci = fits.getdata('f160w_match_drz_sci.fits')\n", - "fig = plt.figure(figsize=(14, 14))\n", - "plt.imshow(sci, vmin=-0.4, vmax=3, cmap='Greys_r', origin='lower')\n", - "plt.colorbar()\n", + "fig = plt.figure(figsize=(10, 10), dpi=130)\n", + "plt.imshow(sci, vmin=0.5, vmax=3, cmap='Greys_r', origin='lower')\n", + "plt.colorbar(shrink=0.85, pad=0.01)\n", "plt.show()" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sci = fits.getdata('f160w_match_drz_wht.fits')\n", + "fig = plt.figure(figsize=(10, 10), dpi=130)\n", + "plt.imshow(sci, vmin=0, vmax=4000, cmap='Greys_r', origin='lower')\n", + "plt.colorbar(shrink=0.85, pad=0.01)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 6. Conclusion \n", + "Thank you for going through this notebook. You should now have all the necessary tools for assessing the
\n", + "appropriate `skymethod` parameter to use when combining images. After completing this notebook you
\n", + "should be more familiar with:
\n", + "    • How to effectively use `astroquery` to download FLT and DRZ files.
\n", + "    • Checking the WCS of images and aligning them to a reference catalog with `TweakReg`.
\n", + "    • Combining data with `AstroDrizzle` taking into account the `skymethod`
\n", + "
\n", + "**Congratulations, you have completed the notebook!**\n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "## Additional Resources \n", + "\n", + "* [DrizzlePac Handbook](https://hst-docs.stsci.edu/drizzpac)\n", + "* [HST Help Desk](https://stsci.service-now.com/hst)\n", + "* [Other DrizzlePac Notebooks](https://github.com/spacetelescope/hst_notebooks/tree/main/notebooks/DrizzlePac)" + ] + }, { "cell_type": "markdown", "metadata": {}, "source": [ + "## About this Notebook \n", "\n", + "Created: 14 Dec 2018; C. Martlin & J. Mack
\n", + "Updated: 16 Nov 2023; K. Huynh & J. Mack
\n", + "Updated: 23 Jul 2024; J. Mack & B. Kuhn
\n", "\n", - " Author: C. Martlin & K.Huynh, STScI WFC3 Team\n", - " Created: December 14, 2018\n", - " Updated: December 22, 2023\n", - " " + "**Source:** https://github.com/spacetelescope/hst_notebooks
" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], - "source": [] + "source": [ + "## Citations \n", + "If you use Python packages for published research, please cite the authors. Follow these links for more
\n", + "information about citing packages such as `astropy`, `astroquery`, `matplotlib`, `pandas`, etc.:
\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/6325653)
\n", + "* [Citing `matplotlib`](https://matplotlib.org/stable/users/project/citing.html)
\n", + "* [Citing `numpy`](https://numpy.org/citing-numpy/)
\n", + "* [Citing `pandas`](https://pandas.pydata.org/about/citing.html)
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[Top of Page](#top)\n", + "\"Space " + ] } ], "metadata": { @@ -571,7 +954,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.7" + "version": "3.12.4" }, "varInspector": { "cols": {