{ "cells": [ { "cell_type": "markdown", "id": "920b5c43", "metadata": {}, "source": [ "\n", "# STIS DrizzlePac Tutorial " ] }, { "cell_type": "markdown", "id": "a77d9e71", "metadata": { "toc": true }, "source": [ "

Table of Contents

\n", "
" ] }, { "cell_type": "markdown", "id": "923edb67", "metadata": {}, "source": [ "_____________________________\n", "## Learning Goals\n", "\n", "By the end of this tutorial, you will:\n", "- Learn how to download STIS data from [MAST](https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html) using [`astroquery`](https://astroquery.readthedocs.io/en/latest/mast/mast.html).\n", "- CTI correct STIS CCD data with the new pixel-based CTI code [`stis_cti`](https://stis-cti.readthedocs.io) (see [webpage](https://www.stsci.edu/hst/instrumentation/stis/data-analysis-and-software-tools/pixel-based-cti) for more details).\n", "- Align images to sub-pixel accuracy for all three STIS detectors (CCD, NUV MAMA, FUV MAMA) using [`tweakreg`](https://drizzlepac.readthedocs.io/en/latest/tweakreg.html) from [DrizzlePac](https://drizzlepac.readthedocs.io/en/latest/index.html).\n", "- Combine images using [`astrodrizzle`](https://drizzlepac.readthedocs.io/en/latest/astrodrizzle.html) from [DrizzlePac](https://drizzlepac.readthedocs.io/en/latest/index.html).\n", "\n", "## Introduction\n", "\n", "The [Space Telescope Imaging Spectrograph (STIS)](https://www.stsci.edu/hst/instrumentation/stis) instrument on board the _Hubble Space Telescope_ (HST) has three detectors: a charge-coupled device (CCD) and two multi-anode microchannel array (MAMA) detectors for near- and far-ultraviolet wavelengths (NUV MAMA and FUV MAMA respectively). The CCD detector suffers from charge transfer inefficiency (CTI) which can be corrected for directly on the CCD images using a pixel-based algorithm (`stis_cti`; based on work by [Anderson & Bedin 2010](https://ui.adsabs.harvard.edu/abs/2010PASP..122.1035A/abstract), see [webpage](https://www.stsci.edu/hst/instrumentation/stis/data-analysis-and-software-tools/pixel-based-cti)).\n", "\n", "In order to maximize the science capabilities of images taken with STIS, they can be accurately aligned and combined (or 'drizzled'). The code used to do this is from the DrizzlePac package which has been extensively tested on Wide-Field Camera 3 (WFC3) and Advanced Camera for Surveys (ACS) data but was previously unsupported for STIS. This tutorial gives alignment and drizzling examples for STIS images of standard star fields taken for each of its three detectors. We refer users to the [DrizzlePac Handbook](https://www.stsci.edu/scientific-community/software/drizzlepac.html) (Gonzaga et al. 2015, Hoffmann et al. 2021) for more details on the drizzling process and DrizzlePac codes.\n", "\n", "### STIS Data Formats\n", "\n", "**CCD data:** For this example, only post-SM4 CCD data taken on primary science amplifier D are used. These are the best calibrated STIS CCD images, the default for science observations, and can be CTI corrected with the new pixel-based code `stis_cti`. For alignment and drizzling, we use STIS CCD images that have been fully calibrated (sx2.fits images: dark- and bias-subtracted, flat-fielded, sky-subtracted, summed individual cosmic-ray (CR) split images with CR-rejection and distortion corrections already applied) that have also been CTI corrected (s2c.fits images).\n", "\n", "**MAMA data:** MAMA detectors do not require CTI corrections or CR-rejection, and STIS MAMA imaging is typically taken in single exposures (i.e., not CR split). For this example, we only use MAMA data with one science extension (i.e., three extensions total: `SCI`, `ERR`, `DQ`) as is typical for STIS MAMA data. If the MAMA data has more than one science extension (`NEXTEND>3`), the science extensions can be summed and saved as a new fits file with this as the first (and only) science extension. These should also be available as sfl.fits images from MAST where this is relevant. For alignment and drizzling, we use STIS MAMA images that have been fully calibrated (x2d.fits: dark- and bias-subtracted, flat-fielded, sky-subtracted, distortion corrections already applied).\n", "\n", "Typically `tweakreg` and `astrodrizzle` require individual dark- and bias-subtracted, flat-fielded images (flt.fits), as for WFC3 and ACS. For STIS data, a few workarounds are required to ensure these codes can run on the STIS data. Additionally, specific parameters are set to ensure that additional calibrations (e.g, sky subtraction, CR-rejection, distortion corrections) are not applied to the already calibrated STIS data.\n", "\n", "See the [STIS Instrument Handbook](https://hst-docs.stsci.edu/stisihb) and [STIS Data Handbook](https://hst-docs.stsci.edu/stisdhb) (specifically [Section 2.2](https://hst-docs.stsci.edu/stisdhb/chapter-2-stis-data-structure/2-2-types-of-stis-files) on file types) for more information.\n", "____________________\n", "## Imports\n", "\n", "- _os_ for operating system commands\n", "- _sys_ to set system preferences\n", "- _glob_ for getting file names\n", "- _shutil_ for copying files\n", "- _subprocess_ for running system commands in the Notebook\n", "- _numpy_ to handle array functions\n", "- _pandas_ for dataframe storage and manipulation\n", "- _matplotlib.pyplot_ for plotting images\n", "- _tempfile TemporaryDirectory_ for creating temporary storage directories while downloading and sorting data\n", "- _astropy.io fits_ for working with FITS files\n", "- _astropy.io ascii_ for working with ascii files\n", "- _astropy.table Table_ for creating tidy tables\n", "- _astropy.visualization ZScaleInterval_ for scaling images for display\n", "- _stis_cti_ for CTI correcting STIS data\n", "- _drizzlepac_ for alignment (_tweakreg_) and image drizzling (_astrodrizzle_)\n", "- _astroquery_ for querying and downloading data from MAST\n", "- _copy_files_ custom function (in directory) for creating directories, checking for data and copying files" ] }, { "cell_type": "code", "execution_count": null, "id": "6f4e95b8", "metadata": {}, "outputs": [], "source": [ "# Load packages\n", "import os\n", "import sys\n", "import glob\n", "import shutil\n", "import subprocess\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib\n", "import matplotlib.pyplot as plt\n", "from tempfile import TemporaryDirectory\n", "from astropy.io import fits\n", "from astropy.io import ascii\n", "from astropy.table import Table\n", "from astropy.visualization import ZScaleInterval\n", "\n", "# STScI packages\n", "from drizzlepac import tweakreg\n", "from drizzlepac import astrodrizzle as ad\n", "from astroquery.mast import Observations\n", "\n", "# Custom package, copy_files.py should be in the same directory\n", "import copy_files as cf\n", "\n", "# STIS packages\n", "from stis_cti import stis_cti, archive_dark_query\n", "from stis_cti import __version__ as stis_cti_version\n", "\n", "# Notebook settings\n", "from IPython.core.interactiveshell import InteractiveShell\n", "InteractiveShell.ast_node_interactivity = \"all\"\n", "\n", "print(f'stis_cti v{stis_cti_version}')" ] }, { "cell_type": "code", "execution_count": null, "id": "e00516cb", "metadata": {}, "outputs": [], "source": [ "# Set plotting defaults\n", "matplotlib.rcParams['figure.figsize'] = [15, 5]\n", "matplotlib.rcParams['image.origin'] = 'lower'\n", "matplotlib.rcParams['image.aspect'] = 'equal'\n", "matplotlib.rcParams['image.cmap'] = 'inferno'\n", "matplotlib.rcParams['image.interpolation'] = 'none'\n", "%matplotlib inline\n", "\n", "# Setting pandas display options\n", "pd.set_option('display.expand_frame_repr', False)\n", "pd.set_option('display.max_rows', 3000)\n", "pd.set_option('display.max_columns', 500)\n", "pd.set_option('display.width', 1000)\n", "\n", "# Set numpy array display options\n", "np.set_printoptions(threshold=sys.maxsize)" ] }, { "cell_type": "markdown", "id": "0917a27a", "metadata": {}, "source": [ "__________________\n", "## Set & Make Directories\n", "Set root directory for project (`root_dir`), a reference file directory (`ref_dir`, explained in more detail below but set this path to any existing reference file directories if present), and a cache directory for temporary data storage for downloads (``cache_dir``). The remaining directories will be nested within the root directory and are created below. " ] }, { "cell_type": "code", "execution_count": null, "id": "562dad1d", "metadata": {}, "outputs": [], "source": [ "def check_mkdir(new_dir):\n", " \"\"\"Check for and make directory.\n", " \n", " Function to check for a directory and make it if it doesn't exist.\n", " \n", " Parameters\n", " ----------\n", " new_dir : str\n", " Path of directory to check for/make.\n", " \n", " \"\"\"\n", " \n", " # Check for directory\n", " if os.path.exists(new_dir):\n", " print('Directory exists: {}'.format(new_dir))\n", " else:\n", " # Make directory\n", " os.makedirs(new_dir, 0o774)\n", " print('Created directory: {}'.format(new_dir))" ] }, { "cell_type": "code", "execution_count": null, "id": "a25c4e92", "metadata": {}, "outputs": [], "source": [ "# Get the Notebook current working directory to handle relative paths\n", "cwd = os.getcwd()\n", "\n", "# Set temporary data download store\n", "cache_dir = './data_cache' # Update this as needed\n", "\n", "# Set the reference file directory\n", "ref_dir = './reference_files' # Update this if another reference file directory exists\n", "\n", "# Set root and data directories\n", "root_dir = './drizpac' # Update this as needed\n", "dat_dir = os.path.join(root_dir, 'data')\n", "ccd_dir = os.path.join(dat_dir, 'ccd_data')\n", "mama_dir = os.path.join(dat_dir, 'mama_data')\n", "\n", "# Set directories for CTI correction of CCD data\n", "cti_dir = os.path.join(dat_dir, 'ccd_cti')\n", "science = os.path.join(cti_dir, 'science')\n", "darks = os.path.join(cti_dir, 'darks')\n", "ref_cti = os.path.join(cti_dir, 'ref')\n", "\n", "# Set directories for alignment with tweakreg\n", "ccd_twk = os.path.join(dat_dir, 'ccd_twk')\n", "nuv_twk = os.path.join(dat_dir, 'nuv_twk')\n", "fuv_twk = os.path.join(dat_dir, 'fuv_twk')\n", "\n", "# Set directories for drizzling with astrodrizzle\n", "ccd_drz = os.path.join(dat_dir, 'ccd_drz')\n", "nuv_drz = os.path.join(dat_dir, 'nuv_drz')\n", "fuv_drz = os.path.join(dat_dir, 'fuv_drz')\n", "\n", "# Check for directories and make them if they don't exist\n", "for d in [cache_dir, ref_dir, root_dir, dat_dir, ccd_dir, mama_dir, cti_dir, science, darks, ref_cti, ccd_twk, fuv_twk, nuv_twk, ccd_drz, fuv_drz, nuv_drz]:\n", " check_mkdir(d)" ] }, { "cell_type": "markdown", "id": "572c258b", "metadata": {}, "source": [ "## Setting Up Reference File Directory\n", "\n", "STScI codes (e.g., `stis_cti`, `DrizzlePac`) need to access reference files to ensure the HST data they run on are properly calibrated. A local store for the reference files is required and a systerm variable should be set for the codes to access those files. The directory variable for STIS data is called \"`oref`\" by convention. \n", "\n", "We create the `oref` environment variable first and populate it with the necessary reference files later in the Notebook. If you already have the required STIS reference files in an existing `oref` directory, set the `oref` environment variable to that directory path.\n", "\n", "We can assign a system variable in three different ways, depending on whether we are working from:\n", "\n", "1. The command line\n", "2. A Python/Jupyter environment\n", "3. A Jupyter Notebook\n", "\n", "|Unix-style Command Line| Python/Jupyter | Jupyter Notebook|\n", "|-|-|-|\n", "| export oref='./reference_files/references/...' | os.environ[\\'oref\\'] = \\'./reference_files/references/...\\' | %env oref ./reference_files/references/... |\n", "\n", "Note that this system variable must be set again with every new instance of a terminal or Notebook or the `export` command can be added to your `.bash_profile` or equivalent file.\n", "\n", "The [Calibration Reference Data System (CRDS)](https://hst-crds.stsci.edu/) both stores the reference files and determines the mapping of reference files to observations. The `crds` tool can find, update, and download the best reference files for a particular observation. The [documentation for `crds`](https://hst-crds.stsci.edu/static/users_guide/index.html) describes many of the more advanced options. In this Notebook we provide a demonstration on how to obtain updated reference file information stored in the FITS headers of observations and download those files to a local reference file directory.\n", "\n", "First the following environment variables are needed:\n", "- `CRDS_SERVER_URL`: Location of the CRDS server.\n", "- `CRDS_PATH`: Path to where reference files will be downloaded.\n", "- `oref`: Path that the reference files will be downloaded to within `CRDS_PATH`" ] }, { "cell_type": "code", "execution_count": null, "id": "959bbbcd", "metadata": {}, "outputs": [], "source": [ "# Set environment variables\n", "os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'\n", "os.environ['CRDS_PATH'] = os.path.abspath(ref_dir)\n", "os.environ['oref'] = os.path.join(os.path.abspath(ref_dir), 'references', 'hst', 'stis') + os.path.sep # Trailing slash important\n", "# print('oref: {}'.format(os.environ['oref']))" ] }, { "cell_type": "markdown", "id": "9e1f2ce3", "metadata": {}, "source": [ "## Downloading Data\n", "The data used in this example are from STIS monitoring observations of two standard star fields: NGC 5139 for the CCD and NGC 6681 for the MAMAs. We use two programs for each field from 2010 and 2011. The data are taken at the same position angle (PA) with some dithering and in the following filters:\n", "- CCD: 50CCD\n", "- NUV MAMA: F25SRF2, F25QTZ, F25CN182\n", "- FUV MAMA: F25SRF2, F25QTZ, 25MAMA\n", "\n", "The data are downloaded from [MAST](https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html) using [`astroquery`](https://astroquery.readthedocs.io/en/latest/mast/mast.html), ~288 MB of data for the CCD and ~437 MB of data for the MAMAs." ] }, { "cell_type": "code", "execution_count": null, "id": "5cb21e15", "metadata": {}, "outputs": [], "source": [ "# Set proposal IDs (PIDs)\n", "ccd_pids = [11854, 12409]\n", "mama_pids = [11856, 12413] # Both NUV & FUV MAMA observations taken in a single PID" ] }, { "cell_type": "code", "execution_count": null, "id": "88edc45c", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Check astroquery Observations options\n", "Observations.get_metadata(\"observations\")\n", "Observations.get_metadata(\"products\")" ] }, { "cell_type": "code", "execution_count": null, "id": "35f126f8", "metadata": {}, "outputs": [], "source": [ "def get_oids(pids):\n", " \"\"\"Get observation IDs from proposal IDs.\n", " \n", " From a list of input proposal IDs (``pids``) get a list of observation IDs \n", " (``oids``) using astroquery.\n", "\n", " Parameters\n", " ----------\n", " pids : list or array_like of str or int\n", " List or array of proposal IDs to find observation IDs for\n", " \n", " Returns\n", " -------\n", " oids : list of str\n", " List of observation IDs within the input proposal IDs\n", " \n", " \"\"\"\n", " \n", " oids = []\n", " \n", " # For each PID get obs IDs\n", " for pid in pids:\n", " obs = Observations.query_criteria(proposal_id='{}'.format(pid)) \n", " products = Observations.get_product_list(obs)\n", " oids += list(np.unique(products['obs_id']))\n", " \n", " print('{} observation IDs found for {} proposal IDs'.format(len(oids), len(pids)))\n", " return oids\n", "\n", "\n", "# Get the obs_id values from the PIDs\n", "ccd_oids = get_oids(ccd_pids)\n", "mama_oids = get_oids(mama_pids)" ] }, { "cell_type": "code", "execution_count": null, "id": "907c5316", "metadata": {}, "outputs": [], "source": [ "def download_data(ids, destination, product_types=None, cache_dir=cache_dir):\n", " '''Downloads MAST data products into a flattened location.\n", " \n", " Downloads data products (``product_type``) from input observation IDs (``ids``) \n", " from MAST and copies them to a single directory (``destination``) from the \n", " temporary download directory (``cache_dir``). Similar to stisteam.getdata(). \n", " Written by Sean Lockwood.\n", " \n", " Parameters\n", " ----------\n", " ids : list of str\n", " List of observation IDs to download data products from\n", " destination : str\n", " Full path to final destination directory for data\n", " product_types : list of str, optional\n", " Names of product types to download for each observation ID (default is None, \n", " means all data products will be downloaded)\n", " cache_dir : str, optional\n", " Full path to temporary data download storage directory (default \n", " ``cache_dir`` as defined above)\n", " \n", " '''\n", " \n", " assert os.path.isdir(destination), 'Destination must be a directory'\n", " \n", " print('\\nDownloading & copying data to {}\\n'.format(destination))\n", "\n", " # Get data products for each observation ID\n", " obs = Observations.query_criteria(obs_id=ids) \n", " products = Observations.get_product_list(obs)\n", " if product_types is not None:\n", " products = products[[x.upper() in product_types for x in products['productSubGroupDescription']]]\n", " \n", " # Download data and combine into the destination directory\n", " with TemporaryDirectory(prefix='downloads', dir=cache_dir) as d:\n", " dl = Observations.download_products(products, mrp_only=False, download_dir=d)\n", " for filename in dl['Local Path']:\n", " shutil.copy(filename, destination)" ] }, { "cell_type": "code", "execution_count": null, "id": "0ac47397", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Download and combine CCD data (~288 MB)\n", "download_data(ccd_oids, ccd_dir, product_types={'RAW', 'EPC', 'SPT', 'ASN', 'SX2'})\n", "\n", "# Download and combine MAMA data (~437 MB)\n", "download_data(mama_oids, mama_dir, product_types={'X2D'})" ] }, { "cell_type": "markdown", "id": "01ab6aa0", "metadata": {}, "source": [ "## Download Reference Files\n", "\n", "The reference files for the STIS images are then updated in the image headers and downloaded (~250 MB). While the `crds.bestrefs` tool is also accessible inside of Python, it was designed with a command line interface in mind, therefore it is easiest to use it this way." ] }, { "cell_type": "code", "execution_count": null, "id": "4fc2d435", "metadata": {}, "outputs": [], "source": [ "# Find and update the headers with the best reference files and download them\n", "def download_ref_files(image_list):\n", " print('Downloading reference files for:\\n{}\\n'.format(image_list))\n", " images = ' '.join(image_list)\n", " subprocess.check_output('crds bestrefs --files {} --sync-references=1 --update-bestrefs'.format(images), shell=True, stderr=subprocess.DEVNULL)\n", "\n", "\n", "# Check in the Notebook directory\n", "os.chdir(cwd)\n", "\n", "# Get a list of images for the CCD and download reference files\n", "ccd_list = glob.glob(os.path.join(ccd_dir, '*_raw.fits'))\n", "download_ref_files(ccd_list)\n", "\n", "# Get a list of images for the MAMAs and download reference files\n", "mama_list = glob.glob(os.path.join(mama_dir, '*_x2d.fits'))\n", "download_ref_files(mama_list)" ] }, { "cell_type": "markdown", "id": "b21443fd", "metadata": {}, "source": [ "## CTI Correct CCD Images\n", "\n", "A pixel-based CTI correction is applied to the STIS CCD data with [`stis_cti`](https://stis-cti.readthedocs.io). At present, this code can only be run on post-Servicing Mission 4 (SM4 in 2009) CCD data taken on the default science amplifier D. For all other CCD data, the previous empricial method ([`stistools.ctestis`](https://stistools.readthedocs.io/en/latest/ctestis.html)) can be used.\n", "\n", "First the necessary files are copied over to the CTI correction directory. Then the darks needed to correct the data are determined and downloaded (about ~1.2 GB of data). Then the `stis_cti` code is run, see this [webpage](https://stis-cti.readthedocs.io/en/latest/stis_cti.html#stis-cti-stis-cti) for a full description of the parameters. The `num_processes` parameter is the maximum number of parallel processes to use when running `stis_cti`, adjust this as needed for your machine." ] }, { "cell_type": "code", "execution_count": null, "id": "cb874b7b", "metadata": {}, "outputs": [], "source": [ "# Set file extensions for CCD image data to copy\n", "exts = ['*_raw.fits', '*_epc.fits', '*_spt.fits', '*_asn.fits']\n", "\n", "# Copy over all CCD files for each extension to the CTI correction directory\n", "os.chdir(cwd)\n", "for ext in exts:\n", " cf.copy_files_check(ccd_dir, science, files='{}'.format(ext))" ] }, { "cell_type": "code", "execution_count": null, "id": "4c8060dd", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Determine dark exposures required for CCD data and download them (~1.2 GB)\n", "needed_darks = archive_dark_query(glob.glob(os.path.join(science, '*_raw.fits')))\n", "\n", "dark_exposures = set()\n", "for anneal in needed_darks:\n", " for dark in anneal['darks']:\n", " dark_exposures.add(dark['exposure'])\n", "\n", "download_data(dark_exposures, darks, product_types={'FLT', 'EPC', 'SPT'})" ] }, { "cell_type": "code", "execution_count": null, "id": "af2339a8", "metadata": {}, "outputs": [], "source": [ "# Run the CTI code on the CCD data\n", "stis_cti(\n", " science_dir=os.path.abspath(science),\n", " dark_dir=os.path.abspath(darks),\n", " ref_dir=os.path.abspath(ref_cti),\n", " num_processes=15,\n", " clean=True,\n", " verbose=2,)" ] }, { "cell_type": "markdown", "id": "5fd56006", "metadata": {}, "source": [ "Run checks on and inspect the CCD CTI corrected images." ] }, { "cell_type": "code", "execution_count": null, "id": "5c894845", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Print information for each file, check those that have been CTI corrected\n", "os.chdir(cwd)\n", "files = glob.glob(os.path.join(os.path.abspath(science), '*raw.fits'))\n", "files.sort()\n", "\n", "psm4 = 0\n", "ampd = 0\n", "ampo = 0\n", "for f in files:\n", " \n", " # Open file header\n", " hdr = fits.getheader(f, 0)\n", " \n", " # Get year obs\n", " yr = int(hdr['TDATEOBS'].split('-')[0])\n", " \n", " # Get file info\n", " if yr <= 2004: \n", " print('\\n...PRE-SM4 ({}), NOT CTE CORRECTED...'.format(yr))\n", " psm4 += 1\n", " elif os.path.exists(f.replace('raw.fits', 'cte.fits')): \n", " print('\\n***CTE CORRECTED, AMP {} ({})***'.format(hdr['CCDAMP'], yr))\n", " ampd += 1\n", " else: \n", " ampo += 1\n", " print('\\n~NON-CTE CORRECTED, AMP {} ({})~'.format(hdr['CCDAMP'], yr))\n", " print('FILE: {}, PID: {}, ROOT: {}'.format(os.path.basename(f), hdr['PROPOSID'], hdr['ROOTNAME']))\n", " print('INST: {}, DETECTOR: {}, AP: {}'.format(hdr['INSTRUME'], hdr['DETECTOR'], hdr['APERTURE']))\n", " print('DATE OBS:{}, PROCESSED: {}'.format(hdr['TDATEOBS'], hdr['DATE']))" ] }, { "cell_type": "code", "execution_count": null, "id": "4a672ef0", "metadata": {}, "outputs": [], "source": [ "# Run checks on data\n", "nraw = len(glob.glob(os.path.join(science, '*raw.fits')))\n", "ncte = len(glob.glob(os.path.join(science, '*cte.fits')))\n", "print('{} raw CCD files, {} CTE corrected files'.format(nraw, ncte)) # Should be equal if all data taken on amp D and post-SM4\n", "print('{} pre-SM4, {} amp D (CTE corr), {} other amp (non-CTE)'.format(psm4, ampd, ampo))" ] }, { "cell_type": "code", "execution_count": null, "id": "83a46227", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Display example CCD images\n", "sci_root = 'obat01050'\n", "sx2_file = os.path.join(ccd_dir, '{}_sx2.fits'.format(sci_root)) # CR rejected, distortion corrected\n", "s2c_file = os.path.join(science, '{}_s2c.fits'.format(sci_root)) # CR rejected, distortion corrected, CTI corrected\n", "\n", "ccd_files = [sx2_file]\n", "cti_files = [s2c_file]\n", "\n", "for ccdf, ctif in zip(ccd_files, cti_files):\n", " \n", " # Plot whole image\n", " fig, axes = plt.subplots(1, 2)\n", " _ = fig.set_size_inches(15, 7)\n", "\n", " vmin, vmax = -10, 10\n", " for filename, ax in zip([ccdf, ctif], axes):\n", " _ = ax.imshow(fits.getdata(filename, ext=1), vmin=vmin, vmax=vmax)\n", " _ = ax.set_title(os.path.basename(filename))\n", " plt.show()\n", "\n", " # Zoom in near the bottom center\n", " fig, axes = plt.subplots(1, 2)\n", " _ = fig.set_size_inches(15, 7)\n", "\n", " vmin, vmax = -10, 10\n", " for filename, ax in zip([ccdf, ctif], axes):\n", " _ = ax.imshow(fits.getdata(filename, ext=1), vmin=vmin, vmax=vmax)\n", " _ = ax.set_title('{} zoom'.format(os.path.basename(filename)))\n", " _ = ax.set_xlim(300, 800)\n", " _ = ax.set_ylim(0, 300)\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "69243a83", "metadata": {}, "source": [ "## File Information\n", "\n", "Information is pulled from file headers and combined into a dataframe. The information is used for file sorting and determining appropriate alignment and drizzle parameters. We select only CCD data taken with primary science amplifier D and MAMA data with only one science extension for this example (see the **Introduction - STIS Data Formats** section for more details)." ] }, { "cell_type": "code", "execution_count": null, "id": "e0a52652", "metadata": {}, "outputs": [], "source": [ "# Make dataframe of image properties\n", "os.chdir(cwd)\n", "det_dirs = [science, mama_dir]\n", "exts = ['*_s2c.fits', '*_x2d.fits']\n", "\n", "# Set variable arrays\n", "fnames = []\n", "det_names = []\n", "rnames = []\n", "pids = []\n", "yrs = []\n", "tdats = []\n", "dats = []\n", "dets = []\n", "aps = []\n", "nexs = []\n", "texps = []\n", "amps = []\n", "gains = []\n", "pas = []\n", "ras = []\n", "decs = []\n", "pt1s = []\n", "pt2s = []\n", "\n", "# Loop over folders and file extensions\n", "for det_dir, ext in zip(det_dirs, exts):\n", " \n", " # Get list of file names\n", " files = glob.glob(os.path.join(det_dir, ext))\n", " files.sort()\n", " \n", " # Loop over each file\n", " for f in files:\n", " # Open file header\n", " hdr = fits.getheader(f, 0) + fits.getheader(f, 1)\n", "\n", " # Get image properties\n", " fname = os.path.basename(f)\n", " rname = hdr['ROOTNAME'].strip()\n", " pid = hdr['PROPOSID']\n", " yr = int(hdr['TDATEOBS'].split('-')[0])\n", " tdat = hdr['TDATEOBS'].strip()\n", " dat = hdr['DATE'].strip().split('T')[0]\n", " det = hdr['DETECTOR'].strip()\n", " nex = hdr['NEXTEND']\n", " texp = hdr['TEXPTIME']\n", " ap = hdr['APERTURE'].strip()\n", " pa = hdr['PA_APER']\n", " ra = hdr['RA_TARG']\n", " dec = hdr['DEC_TARG']\n", " pt1 = hdr['POSTARG1']\n", " pt2 = hdr['POSTARG2']\n", "\n", " # For CCD data, store different header information\n", " if 'CCD' in det:\n", " amp = hdr['CCDAMP'].strip()\n", " gain = hdr['CCDGAIN']\n", " \n", " # For CCD data, only using amp D data (i.e., the only data that is CTI corrected)\n", " if 'D' in amp:\n", " if ext == '*_s2c.fits':\n", " det_name = 'ccd_cti'\n", " else:\n", " det_name = 'ccd'\n", " else: \n", " # Name excluded CCD files\n", " if ext == '*_s2c.fits':\n", " det_name = 'xccd_cti'\n", " else:\n", " det_name = 'xccd'\n", "\n", " # Store detector name for MAMAs\n", " elif 'MAMA' in det: \n", " \n", " # Take only single SCI exposure MAMA files (i.e. 3 extensions, SCI, ERR, DQ)\n", " if nex == 3:\n", " if 'FUV-MAMA' in det:\n", " det_name = 'fuv_mama'\n", " elif 'NUV-MAMA' in det:\n", " det_name = 'nuv_mama'\n", " else: \n", " if 'FUV-MAMA' in det:\n", " det_name = 'xfuv_mama'\n", " elif 'NUV-MAMA' in det:\n", " det_name = 'xnuv_mama'\n", " \n", " # MAMAs do not have amp or gain info\n", " amp = '-'\n", " gain = '-'\n", " \n", " # Save file info\n", " fnames.append(fname)\n", " det_names.append(det_name)\n", " rnames.append(rname)\n", " pids.append(pid)\n", " yrs.append(yr)\n", " tdats.append(tdat)\n", " dats.append(dat)\n", " dets.append(det)\n", " aps.append(ap)\n", " nexs.append(nex)\n", " texps.append(texp)\n", " amps.append(amp)\n", " gains.append(gain)\n", " pas.append(pa)\n", " ras.append(ra)\n", " decs.append(dec)\n", " pt1s.append(pt1)\n", " pt2s.append(pt2)\n", "\n", "# Create data frame of image properties\n", "df = pd.DataFrame({})\n", "df['file'] = fnames\n", "df['type'] = det_names\n", "df['rootname'] = rnames\n", "df['pid'] = pids\n", "df['year_obs'] = yrs\n", "df['date_obs'] = tdats\n", "df['date_proc'] = dats\n", "df['detector'] = dets\n", "df['aperture'] = aps\n", "df['nextend'] = nexs\n", "df['texptime'] = texps\n", "df['amp'] = amps\n", "df['gain'] = gains\n", "df['pa_aper'] = pas\n", "df['ra'] = ras\n", "df['dec'] = decs\n", "df['postarg1'] = pt1s\n", "df['postarg2'] = pt2s\n", "df = df.sort_values(by=[\"type\", \"date_obs\"]).reset_index(drop=True)\n", "df" ] }, { "cell_type": "code", "execution_count": null, "id": "242f4012", "metadata": {}, "outputs": [], "source": [ "# Get exposure times to aid with determining alignment parameters\n", "print('\\nCCD-CTI CORRECTED EXP TIMES ({} files)'.format(len(df.loc[df['type'] == 'ccd_cti'])))\n", "np.sort(df['texptime'].loc[df['type'] == 'ccd_cti'].unique())\n", "\n", "print('\\nNUV MAMA EXP TIMES ({} files)'.format(len(df.loc[df['type'] == 'nuv_mama'])))\n", "np.sort(df['texptime'].loc[df['type'] == 'nuv_mama'].unique())\n", "\n", "print('\\nFUV MAMA EXP TIMES ({} files)'.format(len(df.loc[df['type'] == 'fuv_mama'])))\n", "np.sort(df['texptime'].loc[df['type'] == 'fuv_mama'].unique())" ] }, { "cell_type": "code", "execution_count": null, "id": "0014e069", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Check for directories, make them if they don't exist, copy files to alignment directories\n", "os.chdir(cwd)\n", "cf.copy_files_check(science, ccd_twk, files=df['file'].loc[df['type'] == 'ccd_cti']) # s2c: Distortion+CR+CTI corrected CCD\n", "cf.copy_files_check(mama_dir, nuv_twk, files=df['file'].loc[df['type'] == 'nuv_mama']) # x2d: Distortion corrected NUV-MAMA\n", "cf.copy_files_check(mama_dir, fuv_twk, files=df['file'].loc[df['type'] == 'fuv_mama']) # x2d: Distortion corrected FUV-MAMA" ] }, { "cell_type": "markdown", "id": "0bcfb259", "metadata": {}, "source": [ "## Image Alignment\n", "Use `tweakreg` to align the images from the three STIS detectors to sub-pixel accuracy.\n", "\n", "### Tips and Tricks for Using [tweakreg](https://drizzlepac.readthedocs.io/en/latest/tweakreg.html)\n", "\n", "- For `tweakreg` to work on the already distortion-corrected STIS data, it needs the missing header keyword `IDCSCALE` ('default' plate scale for the detector) as done below.\n", "- `tweakreg` can align images to better than 0.1 pixel accuracy but it can vary between ~0.05-0.5 pixels depending on the image/number of sources/relative alignment to the reference image etc.\n", "- **NOTE:** `tweakreg` changes the image headers after each run, unless `updatehdr=False` is set. Therefore, for testing alignment parameters, always use `updatehdr=False` until happy with the parameters. For the final run, set `updatehdr=True` to update the image headers with their new WCS solution ready for drizzling.\n", "- Select a good quality reference image, e.g., a deeper exposure, roughly centered within the area covered by the other images, no obvious issues/contamination.\n", "- Testing of parameters should, to first order, improve the accuracy of alignment (`XRMS/YRMS` parameters). But don't rely on this value alone to verify the quality of fit.\n", "- Pixel shifts (`XSH/YSH` parameters) between images can range between 0-30 pixels (maybe up to 50 pixels). Expected shifts can be sanity checked by eye with an image viewer (e.g. DS9). Very large shifts (100s pixels) likely indicate a false solution.\n", "- The number of sources found in images and used for alignment is less important than the quality of the sources used. See the **Testing Alignment Parameters** section below for code to overplot the sources used for alignment on the images.\n", "- Below are common parameters to adjust to align images ([input image params](https://drizzlepac.readthedocs.io/en/latest/imagefindpars.html)/[reference image params](https://drizzlepac.readthedocs.io/en/latest/refimagefindpars.html)), trial and error is best for honing in on a solution:\n", " - `threshold`: sigma limit above background for source detection.\n", " - Note that sigma is highly image dependent and can vary widely (e.g. from 0.5 to 200 for the same image and get similar results). It is worth exploring a wide parameter space to find an initial solution.\n", " - `minobj`: Minimum number of objects to be matched between images. \n", " - 10 is a reasonable lower limit, sometimes lower is needed (e.g. for FUV MAMA with fewer sources) but be warned of going too low and false detections/solutions will be found\n", " - `peakmax`: Sets a maximum value of sources in the image. \n", " - Useful for selecting bright stars (needed in FUV-MAMA with so few sources) or avoiding saturated stars.\n", " - `use_sharp_round=True`: Select well-defined sources.\n", " - `conv_width`: Convultion kernel width, recommended 2x the PSF FWHM. \n", " - Sometimes for lower quality/blurrier images this has to be increased to find a solution.\n", " - `searchrad`: Search radius for a match (default units arcsec). \n", " - Sometimes you need to increase this by quite a bit to find a solution. Default 1, but you can use up to 10-20 depending on images.\n", "- There are many more parameters to explore but reasonable alignment solutions should be obtainable for most STIS images with some combination of these.\n", "- If you need to re-run `tweakreg` on already tweaked images (ran with `updatehdr=True`), make a fresh directory and copy over clean files again using the cell above (`copy_files_check`)." ] }, { "cell_type": "code", "execution_count": null, "id": "5bdff76a", "metadata": {}, "outputs": [], "source": [ "# Add in missing header keyword 'IDCSCALE' for tweakreg to work on STIS images\n", "\n", "# Set tweak directories, file extensions and detector names\n", "twk_dirs = [ccd_twk, nuv_twk, fuv_twk]\n", "exts = ['*_s2c.fits', '*_x2d.fits', '*_x2d.fits']\n", "dets = ['CCD', 'NUV MAMA', 'FUV MAMA']\n", "\n", "# Loop over each detector\n", "for twk_dir, ext, det in zip(twk_dirs, exts, dets):\n", " os.chdir(cwd)\n", " os.chdir(os.path.abspath(twk_dir))\n", " print('Updating headers for {} files in {}'.format(ext, twk_dir))\n", " \n", " # Loop over each file\n", " for i, f in enumerate(glob.glob(ext)):\n", " hdu = fits.open(f, 'update')\n", " \n", " # Set 'default' plate scale for the detector\n", " if det == 'CCD':\n", " idc = 0.05072\n", " elif det == 'NUV MAMA':\n", " idc = 0.0246037\n", " elif det == 'FUV MAMA':\n", " idc = 0.024395\n", " \n", " # Create keyword in header and save\n", " hdu[1].header['IDCSCALE'] = idc\n", " hdu.close()" ] }, { "cell_type": "markdown", "id": "526ce850", "metadata": {}, "source": [ "### Select Reference Images\n", "Based on the dataframe of image properties in the **File Summary** section, identify an appropriate reference image (`ref_img`) for each detector (see **Tips and Tricks for Using tweakreg** section above). Adjust `tweakreg` parameters per detector to accurately align images onto the reference image. **Set `update=False` for testing and change to `update=True` for final run to update image headers.**" ] }, { "cell_type": "code", "execution_count": null, "id": "9adfeb8d", "metadata": {}, "outputs": [], "source": [ "# Set the reference files for each detector\n", "ccd_ref = 'obat01050_s2c.fits'\n", "nuv_ref = 'obav01v9q_x2d.fits'\n", "fuv_ref = 'obav01w4q_x2d.fits'" ] }, { "cell_type": "code", "execution_count": null, "id": "cd449d0c", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Set the reference image paths\n", "os.chdir(cwd)\n", "cref = os.path.join(ccd_twk, ccd_ref)\n", "nref = os.path.join(nuv_twk, nuv_ref)\n", "fref = os.path.join(fuv_twk, fuv_ref)\n", "\n", "# Plot reference images\n", "fig, axes = plt.subplots(3, 1)\n", "_ = fig.set_size_inches(7, 23)\n", "\n", "vmin, vmax = -10, 10\n", "for filename, ax in zip([cref, nref, fref], axes):\n", " _ = ax.imshow(fits.getdata(filename, ext=1), vmin=vmin, vmax=vmax)\n", " _ = ax.set_title(os.path.basename(filename))\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "c1420f97", "metadata": {}, "source": [ "### CCD (CTI Corrected) Image Alignment\n", "\n", "Parameters tested for a couple of different exposure times are given below: <=10s and 60s. \n", "\n", "Alignment accuracy for these images ranges between ~0.04--0.2 pixels, with average ~0.11 pixels." ] }, { "cell_type": "code", "execution_count": null, "id": "cfcdef60", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# -------------------\n", "# CCD - CTI CORRECTED\n", "# -------------------\n", "\n", "# Keep flag as FALSE until happy with parameters\n", "# Set to TRUE for final run\n", "update = False\n", "\n", "# Move into tweak directory\n", "os.chdir(cwd)\n", "os.chdir(ccd_twk)\n", "print(ccd_twk)\n", "\n", "# Set reference image and extension name\n", "ref_img = ccd_ref\n", "ext = '*_s2c.fits'\n", "\n", "# Set files to align to reference image (remove ref image from list)\n", "files = glob.glob(ext)\n", "files.remove(ref_img)\n", "files.sort()\n", "print('\\nAligning {} {} files to {}: {}\\n'.format(len(files), ext, ref_img, files))\n", "\n", "# Loop over files, get image info, set parameters, align\n", "for i, f in enumerate(files):\n", " \n", " # Get file info\n", " print('\\n+++++++++++++++++++++++++++++++++++++++++++++++++++')\n", " hdr = fits.getheader(f, 0)\n", " texp = hdr['TEXPTIME']\n", " print('FILE ({}/{}): {}, ROOT: {}, TEXP:{}'.format(i+1, len(files), f, hdr['ROOTNAME'], texp))\n", " \n", " # Set params based on exposure time, print details\n", " if texp <= 10:\n", " minobj, ithresh, rthresh, peak, convw = 5, 1.7, 1.0, 1500, 3.5 # tested for <=10s CCD exposures\n", " elif texp == 60:\n", " minobj, ithresh, rthresh, peak, convw = 10, 1.7, 1.0, 900, 3.5 # tested for 60s CCD exposures\n", " print('MINOBJ: {}, THRESH-IMG: {}, THRESH-REF: {}, PEAKMAX:{}, CONVW: {}'.format(minobj, ithresh, rthresh, peak, convw))\n", " print('+++++++++++++++++++++++++++++++++++++++++++++++++++\\n')\n", " \n", " # Run tweakreg for each file\n", " tweakreg.TweakReg(f, updatehdr=update, conv_width=convw, shiftfile=True, outshifts='{}_shifts.txt'.format(hdr['ROOTNAME']),\n", " refimage=ref_img, clean=False, interactive=False, minobj=minobj, \n", " searchrad=15.0, imagefindcfg={'threshold': ithresh, 'peakmax': peak, 'use_sharp_round': True},\n", " refimagefindcfg={'threshold': rthresh, 'peakmax': peak, 'use_sharp_round': True})" ] }, { "cell_type": "markdown", "id": "e6719a31", "metadata": {}, "source": [ "### NUV MAMA Image Alignment\n", "\n", "Parameters tested for a couple of different exposure times are given below as a guide: 300 to <1000s (used here) and >=1000s. For other 300s images spanning 25 years of STIS data (only two programs used here), adjusting the convolution width (e.g., `conv_width=4.5, 5.5, 6.5`) for specific files enabled an alignment solution to be found. One file (obmi01xqq_x2d.fits in `ofiles`) has an offset from the reference image and required different parameters to find an alignment solution.\n", "\n", "Alignment accuracy for these images ranges between ~0.1--0.5 pixels, with average ~0.3 pixels. These values are higher than the CCD partly due to the larger number of sources and that the MAMA data span three filters. " ] }, { "cell_type": "code", "execution_count": null, "id": "c26a8285", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# -------------------\n", "# NUV MAMA\n", "# -------------------\n", "\n", "# Keep flag as FALSE until happy with parameters\n", "# Set to TRUE for final run\n", "update = False\n", "\n", "# Move into tweak directory\n", "os.chdir(cwd)\n", "os.chdir(nuv_twk)\n", "print(nuv_twk)\n", "\n", "# Set reference image and extension name\n", "ref_img = nuv_ref\n", "ext = '*_x2d.fits'\n", "\n", "# Set files to align to reference image (remove ref image from list)\n", "files = glob.glob(ext)\n", "files.remove(ref_img)\n", "files.sort()\n", "ofiles = ['obmi01xqq_x2d.fits'] # file requiring different parameters\n", "print('\\nAligning {} {} files to {}: {}\\n'.format(len(files), ext, ref_img, files))\n", "\n", "# Loop over files, get image info, set parameters, align\n", "for i, f in enumerate(files):\n", " \n", " # Get file info\n", " print('\\n++++++++++++++++++++++++++++++++++++++++++++++++++++++')\n", " hdr = fits.getheader(f, 0)\n", " texp = hdr['TEXPTIME']\n", " print('FILE ({}/{}): {}, ROOT: {}, TEXP:{}'.format(i+1, len(files), f, hdr['ROOTNAME'], texp))\n", " \n", " # Set params based on exposure time, print details\n", " if f in ofiles:\n", " minobj, ithresh, rthresh, peak, convw, sr = 10, 300, 200, 300, 4.5, 20\n", " else:\n", " minobj, ithresh, rthresh, peak, convw, sr = 10, 0.7, 0.5, 500, 5.5, 20 # tested for 300 to <1000s NUV MAMA exposures\n", " # minobj, ithresh, rthresh, peak, convw, sr = 10, 300, 200, 300, 4.5, 10 # tested for >=1000s NUV MAMA exposures\n", " print('MINOBJ: {}, THRESH-IMG: {}, THRESH-REF: {}, PEAKMAX:{}, CONVW: {}'.format(minobj, ithresh, rthresh, peak, convw))\n", " print('++++++++++++++++++++++++++++++++++++++++++++++++++++++\\n')\n", " \n", " # Run tweakreg for each file\n", " tweakreg.TweakReg(f, updatehdr=update, conv_width=convw, shiftfile=True, outshifts='{}_shifts.txt'.format(hdr['ROOTNAME']),\n", " refimage=ref_img, clean=False, interactive=False, minobj=minobj, \n", " searchrad=sr, imagefindcfg={'threshold': ithresh, 'peakmax': peak, 'use_sharp_round': True},\n", " refimagefindcfg={'threshold': rthresh, 'peakmax': peak, 'use_sharp_round': True})" ] }, { "cell_type": "markdown", "id": "6d2caea5", "metadata": {}, "source": [ "### FUV MAMA Image Alignment\n", "\n", "Parameters tested for a couple of different exposure times are given below as a guide: 400s to <1000s (used here), and >=1000s. For other images spanning 25 years of STIS data (only two programs shown here), adjusting the convolution width (e.g., `conv_width=4.5, 5.5, 6.5`) for specific files enabled an alignment solution to be found. Two files (obmi01y8q_x2d.fits, obmi01yaq_x2d.fits in `ofiles`) have an offset from the reference image and required different parameters to find an alignment solution.\n", "\n", "Alignment accuracy for these images ranges between ~0.1--0.4 pixels, with average ~0.18 pixels. These values are slightly higher than the CCD as the MAMA data span three filters and the FUV MAMA has a broader PSF. " ] }, { "cell_type": "code", "execution_count": null, "id": "e287f7ce", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# -------------------\n", "# FUV MAMA\n", "# -------------------\n", "\n", "# Keep flag as FALSE until happy with parameters\n", "# Set to TRUE for final run\n", "update = False\n", "\n", "# Move into tweak directory\n", "os.chdir(cwd)\n", "os.chdir(fuv_twk)\n", "print(fuv_twk)\n", "\n", "# UPDATE: Set reference image and extension names\n", "ref_img = fuv_ref\n", "ext = '*_x2d.fits'\n", "\n", "# Set files to align to reference image (remove ref image from list)\n", "files = glob.glob(ext)\n", "files.remove(ref_img)\n", "files.sort()\n", "ofiles = ['obmi01y8q_x2d.fits', 'obmi01yaq_x2d.fits'] # files requiring different parameters\n", "print('\\nAligning {} {} files to {}: {}\\n\\n'.format(len(files), ext, ref_img, files))\n", "\n", "for i, f in enumerate(files):\n", " \n", " # Get file info\n", " print('\\n\\n\\n++++++++++++++++++++++++++++++++++++++++++++++++++++++')\n", " hdr = fits.getheader(f, 0)\n", " texp = hdr['TEXPTIME']\n", " print('FILE ({}/{}): {}, ROOT: {}, TEXP:{}'.format(i+1, len(files), f, hdr['ROOTNAME'], texp))\n", " \n", " # Set params based on exposure time, print details\n", " if f in ofiles:\n", " minobj, ithresh, rthresh, peak, convw = 5, 100, 100, 3000, 5.5\n", " else:\n", " minobj, ithresh, rthresh, peak, convw = 5, 200, 50, 3000, 5.5 # tested for 400s to <1000s exposures\n", " # minobj, ithresh, rthresh, peak, convw = 10, 200, 50, 8000, 5.5 # tested for >=1000s exposures\n", " print('MINOBJ: {}, THRESH-IMG: {}, THRESH-REF: {}, PEAKMAX:{}, CONVW: {}'.format(minobj, ithresh, rthresh, peak, convw))\n", " print('++++++++++++++++++++++++++++++++++++++++++++++++++++++\\n')\n", " \n", " # Run tweakreg for each file\n", " tweakreg.TweakReg(f, updatehdr=update, conv_width=convw, shiftfile=True, outshifts='{}_shifts.txt'.format(hdr['ROOTNAME']),\n", " refimage=ref_img, clean=False, interactive=False, minobj=minobj, \n", " searchrad=20.0, imagefindcfg={'threshold': ithresh, 'peakmax': peak, 'use_sharp_round': True},\n", " refimagefindcfg={'threshold': rthresh, 'peakmax': peak, 'use_sharp_round': True})" ] }, { "cell_type": "markdown", "id": "48d57c01", "metadata": {}, "source": [ "### Testing Alignment Parameters\n", "Below is useful code for testing `tweakreg` alignment parameters. If an image has been able to be matched (even if `updatehdr=False`) a `_shifts.txt` file will be created. This code can be used to plot the sources used for alignment on the input image. Use this to verify the quality of sources used for alignment." ] }, { "cell_type": "code", "execution_count": null, "id": "0517a96a", "metadata": {}, "outputs": [], "source": [ "# Overplot sources used for alignment (written by Marc Rafelski)\n", "\n", "# Set input directory, reference image and extension names (example for the CCD)\n", "indir = ccd_twk\n", "ref_img = ccd_ref\n", "ext = '*_s2c.fits'\n", "\n", "# Move into tweak directory\n", "os.chdir(cwd)\n", "os.chdir(indir)\n", "print(indir)\n", "\n", "# Set files to align to reference image (remove ref image from list)\n", "files = glob.glob(ext)\n", "files.remove(ref_img)\n", "files.sort()\n", "\n", "# Get rootnames of files\n", "roots = [x.split('_')[0] for x in files]\n", "print(roots)\n", "\n", "# Loop over rootnames\n", "for root in roots:\n", " # Set input file\n", " infile = '{}_shifts.txt'.format(root)\n", "\n", " # Read in table from tweakreg\n", " shift_tab = Table.read(infile, format='ascii.no_header', names=['file', 'dx', 'dy', 'rot', 'scale', 'xrms', 'yrms'])\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.pprint(max_lines=-1)\n", " \n", " # Overplot images with their alignment sources\n", " for image in shift_tab['file']:\n", " _ = plt.figure(figsize=(10, 10))\n", " data = fits.open(image)['SCI', 1].data\n", " zscale = ZScaleInterval()\n", " z1, z2 = zscale.get_limits(data)\n", " _ = plt.imshow(data, cmap='Greys', origin='lower', vmin=z1, vmax=z2)\n", " match_tab = ascii.read(image[0:13]+'_catalog_fit.match')\n", " x_cord, y_cord = match_tab['col11'], match_tab['col12']\n", " _ = plt.scatter(x_cord, y_cord, s=80, edgecolor='r', facecolor='None', label='{} matched sources ({})'.format(len(match_tab), image))\n", " _ = plt.legend(loc='upper right', fontsize=14, framealpha=0.5)\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "45bb1803", "metadata": {}, "source": [ "## Drizzling Images\n", "\n", "Next, the images for each detector are combined using `astrodrizzle`. The high-level concept behind drizzling images is described in detail in [Section 3.2 of the DrizzlePac Handbook](https://hst-docs.stsci.edu/drizzpac/files/60245881/109774911/1/1642098151920/DrizzlePac_Handbook_v2.pdf). \n", "\n", "Setting the appropriate `final_scale` and `final_pixfrac` parameters for your images takes some thought and testing to avoid gaps in the data. The figure below shows a basic example of the native pixel scale (red squares), shrink factor `final_pixfrac` (blue squares) and output final pixel scale `final_scale` (grid on right) in a drizzle. For more details on the `astrodrizzle` input parameters, see the the [DrizzlePac code webpages](https://drizzlepac.readthedocs.io/en/latest/astrodrizzle.html).\n", "\"Drizzle\"\n", "\n", "`astrodrizzle` can be used to increase the pixel sampling of data if images have dithering and different position angles (PAs). For the STIS data used here, all images are at the same PA and therefore sub-sampling the data is not possible. The example for the STIS data shown below adopts the native pixel scale of each detector as the `final_scale` and no fractional scaling down of each pixel (`final_pixfrac=1.0`) prior to drizzling. Drizzle in this context is a useful tool for creating mosaics of images aligned with `tweakreg`.\n", "\n", "To ensure that `astrodrizzle` does not apply additional calibrations to the already calibrated STIS data, Steps 1-6 can be 'switched off' using the appropriate keywords (see [`astrodrizzle` webpage](https://drizzlepac.readthedocs.io/en/latest/astrodrizzle.html) for keyword information). The CCD data are already CR-rejected and the MAMA data don't require this step, therefore many of those first steps (and their associated keywords) are to do with appropriately applying CR-rejection. Therefore the 'minmed' warning is not important for the STIS images.\n", "\n", "Only data in a single filter for each detector are combined in this example: 50CCD (CCD), F25SRF2 (NUV MAMA), 25MAMA (FUV MAMA). STIS data have cleaned edges after the distortion correction is applied in the `calstis` pipeline. These are adaptively cropped out with the `crop_edges` function defined below for each image to avoid gaps in the final mosaics.\n", "\n", "`astrodrizzle` changes the input image files, so it's advisable to first copy the data to a clean directory prior to drizzling (as done below). If you need to repeat the drizzling process, it's good practice to make another clean directory prior to running `astrodrizzle`." ] }, { "cell_type": "code", "execution_count": null, "id": "b7aa6697", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Copy over aligned images to a drizzle directory (creates direc. and copies files), single filters for MAMAs\n", "os.chdir(cwd)\n", "cf.copy_files_check(ccd_twk, ccd_drz, files=df['file'].loc[df['type'] == 'ccd_cti'])\n", "cf.copy_files_check(nuv_twk, nuv_drz, files=df['file'].loc[(df['type'] == 'nuv_mama') & (df['aperture'] == 'F25SRF2')]) # x2d: Distortion corrected NUV-MAMA\n", "cf.copy_files_check(fuv_twk, fuv_drz, files=df['file'].loc[(df['type'] == 'fuv_mama') & (df['aperture'] == '25MAMA')]) # x2d: Distortion corrected FUV-MAMA" ] }, { "cell_type": "code", "execution_count": null, "id": "c92d77c4", "metadata": {}, "outputs": [], "source": [ "def crop_edges(dq, cut=0.2):\n", " \"\"\"Crop edges of distortion-corrected STIS data.\n", " \n", " Function to crop edges of a distortion-corrected STIS image using its data\n", " quality (DQ) array to avoid gaps in the drizzled mosaics. Function makes cuts \n", " at intervals across the array (default fractions of ``cut=0.2``) and determines \n", " the first and last instance of \"good\" data (DQ=0). The good data footprints in \n", " DQ arrays are not perfectly square due to distortion corrections, hence multiple \n", " cuts are made to find the smallest region of good data in both the x and y \n", " directions to avoid partial data rows/columns. Function returns the indices in y \n", " (``ymin``, ``ymax``) and x (``xmin``, ``xmax``) to crop the STIS image data with.\n", " \n", " Parameters\n", " ----------\n", " dq : ndarray\n", " DQ array of image to crop\n", " cut : float, optional\n", " Value between 0 and 1, fractional step sizes across array (default 0.2 means \n", " cuts through data will be at fractions 0.2, 0.4, 0.6, 0.8 of each axis)\n", "\n", " Returns\n", " -------\n", " ymin : int\n", " Lower index of y axis to crop by\n", " ymax : int\n", " Upper index of y axis to crop by\n", " xmin : int\n", " Lower index of x axis to crop by\n", " xmax : int\n", " Upper index of x axis to crop by\n", " \n", " \"\"\"\n", " \n", " # Set storage arrays\n", " ylos = []\n", " yhis = []\n", " xlos = []\n", " xhis = []\n", " \n", " # Make several cuts through the data\n", " for i in np.arange(0, 1, cut)[1:]:\n", " \n", " # Find the index of the row/column\n", " val = int(dq.shape[0]*i)\n", " \n", " # Get start and end indices of the \"good\" data (DQ=0) in x & y \n", " ylos.append(np.where(dq[:, val] == 0)[0][0])\n", " yhis.append(np.where(dq[:, val] == 0)[0][-1])\n", " xlos.append(np.where(dq[val, :] == 0)[0][0])\n", " xhis.append(np.where(dq[val, :] == 0)[0][-1])\n", "\n", " # Get crop indices: smallest region of good data to avoid partial data rows/columns\n", " ymin, ymax, xmin, xmax = np.max(ylos), np.min(yhis)+1, np.max(xlos), np.min(xhis)+1\n", " print('Cropping indices y= {}:{}, x= {}:{}\\n'.format(ymin, ymax, xmin, xmax)) \n", " return ymin, ymax, xmin, xmax" ] }, { "cell_type": "markdown", "id": "3d4f307c", "metadata": {}, "source": [ "### CCD (CTI Corrected) Image Drizzling\n", "\n", "The CCD images are all observed at a common position angle and RA/Dec with small dithers (hence poorer quality edges in the mosaic below). Sub-sampling the images with `astrodrizzle` is not advisable for these programs as reducing the pixel size results in gaps in the data.\n", "\n", "Figure shows the CCD drizzle of NGC 5139 (left) and the individual CCD reference image used for alignment (right).\n", "\"CCD" ] }, { "cell_type": "code", "execution_count": null, "id": "a5a287ed", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Drizzling images together\n", "# Move into drizzle directory\n", "os.chdir(cwd)\n", "os.chdir(ccd_drz)\n", "print(ccd_drz)\n", "\n", "# Set image extension names\n", "ext = '*_s2c.fits'\n", "\n", "# Set files to drizzle\n", "files = glob.glob(ext)\n", "\n", "# Get pixel scales from images and crop data\n", "ps = []\n", "for i, f in enumerate(files):\n", " # Read in HDU list and header\n", " print('\\n{}'.format(f))\n", " hdu = fits.open(f, mode='update')\n", " hdr = hdu[0].header + hdu[1].header\n", " \n", " # Get all image pixel scales (can change slightly between headers)\n", " ps.append(hdr['PLATESC'])\n", " \n", " # Find cropping indices from DQ array\n", " dq = hdu[3].data\n", " ymin, ymax, xmin, xmax = crop_edges(dq)\n", " \n", " # Crop data down in each extension to remove edges with no data (from STIS distortion corrections)\n", " for j in np.arange(3):\n", " hdu[j+1].data = hdu[j+1].data[ymin:ymax, xmin:xmax]\n", " hdu[j+1].header['NAXIS1'] = hdu[j+1].data.shape[1]\n", " hdu[j+1].header['NAXIS2'] = hdu[j+1].data.shape[0]\n", " \n", " # Update primary header array size information\n", " hdu[0].header['SIZAXIS1'] = hdu[j+1].data.shape[1]\n", " hdu[0].header['SIZAXIS2'] = hdu[j+1].data.shape[0]\n", " \n", " # Save changes and close\n", " hdu.close()\n", "\n", "# Set drizzle parameters\n", "fs = np.min(np.array(ps)) # Final pixel scale (arcsec) of the output image, native pixel scale: ~0.050777\n", "fp = 1.0 # Fraction by which to shrink the input pixels prior to drizzling onto output grid\n", "\n", "# Drizzle images together\n", "ad.AstroDrizzle(files, static=False, skysub=False, driz_separate=False, median=False, blot=False, \n", " driz_cr=False, driz_combine=True, clean=True, build=False, preserve=False, \n", " final_scale=fs, final_pixfrac=fp, final_wht_type='ERR', output='ccd_drz')" ] }, { "cell_type": "markdown", "id": "abac322c", "metadata": {}, "source": [ "### NUV MAMA Image Drizzling\n", "\n", "The NUV images are all observed at a common position angle with large dithers (hence different image depths in the mosaic below). Sub-sampling the images with `astrodrizzle` is not advisable for these programs as reducing the pixel size results in gaps in the data.\n", "\n", "Figure shows the NUV drizzle of NGC 6681 (left) and the individual NUV reference image used for alignment (right).\n", "\n", "\"NUV" ] }, { "cell_type": "code", "execution_count": null, "id": "4c78c1df", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Drizzling images together\n", "# Move into drizzle directory\n", "os.chdir(cwd)\n", "os.chdir(nuv_drz)\n", "print(nuv_drz)\n", "\n", "# Set image extension names\n", "ext = '*_x2d.fits'\n", "\n", "# Set files to drizzle\n", "files = glob.glob(ext)\n", "\n", "# Get pixel scales from images and crop data\n", "ps = []\n", "for i, f in enumerate(files):\n", " # Read in HDU list and header\n", " print('{}'.format(f))\n", " hdu = fits.open(f, mode='update')\n", " hdr = hdu[0].header + hdu[1].header\n", " \n", " # Get all image pixel scales (can change slightly between headers)\n", " ps.append(hdr['PLATESC'])\n", " \n", " # Find cropping indices from DQ array\n", " dq = hdu[3].data\n", " ymin, ymax, xmin, xmax = crop_edges(dq)\n", " \n", " # Crop data down in each extension to remove edges with no data (from STIS distortion corrections)\n", " for j in np.arange(3):\n", " hdu[j+1].data = hdu[j+1].data[ymin:ymax, xmin:xmax]\n", " hdu[j+1].header['NAXIS1'] = hdu[j+1].data.shape[1]\n", " hdu[j+1].header['NAXIS2'] = hdu[j+1].data.shape[0]\n", " \n", " # Update primary header array size information\n", " hdu[0].header['SIZAXIS1'] = hdu[j+1].data.shape[1]\n", " hdu[0].header['SIZAXIS2'] = hdu[j+1].data.shape[0]\n", " \n", " # Save changes and close\n", " hdu.close()\n", " \n", "# Set drizzle parameters\n", "fs = np.min(np.array(ps)) # Final pixel scale (arcsec) of the output image, native pixel scale: ~0.02475\n", "fp = 1.0 # Fraction by which to shrink the input pixels prior to drizzling onto output grid\n", "\n", "# Drizzle images together\n", "ad.AstroDrizzle(files, static=False, skysub=False, driz_separate=False, median=False, blot=False, \n", " driz_cr=False, driz_combine=True, clean=True, build=False, preserve=False, \n", " final_scale=fs, final_pixfrac=fp, final_wht_type='ERR', output='nuv_drz')" ] }, { "cell_type": "markdown", "id": "96967f4b", "metadata": {}, "source": [ "### FUV MAMA Image Drizzling\n", "\n", "The FUV images are all observed at a common position angle with large dithers (hence different image depths in the mosaic below). Sub-sampling the images with `astrodrizzle` is not advisable for these programs as reducing the pixel size results in gaps in the data.\n", "\n", "Figure shows the FUV drizzle of NGC 6681 (left) and the individual FUV reference image used for alignment (right).\n", "\"FUV" ] }, { "cell_type": "code", "execution_count": null, "id": "0ffe96cd", "metadata": {}, "outputs": [], "source": [ "# Drizzling images together\n", "# Move into drizzle directory\n", "os.chdir(cwd)\n", "os.chdir(fuv_drz)\n", "print(fuv_drz)\n", "\n", "# Set image extension names\n", "ext = '*_x2d.fits'\n", "\n", "# Set files to drizzle\n", "files = glob.glob(ext)\n", "\n", "# Get pixel scales from images and crop data\n", "ps = []\n", "for i, f in enumerate(files):\n", " # Read in HDU list and header\n", " print('{}'.format(f))\n", " hdu = fits.open(f, mode='update')\n", " hdr = hdu[0].header + hdu[1].header\n", " \n", " # Get all image pixel scales (can change slightly between headers)\n", " ps.append(hdr['PLATESC'])\n", " \n", " # Find cropping indices from DQ array\n", " dq = hdu[3].data\n", " ymin, ymax, xmin, xmax = crop_edges(dq)\n", " \n", " # Crop data down in each extension to remove edges with no data (from STIS distortion corrections)\n", " for j in np.arange(3):\n", " hdu[j+1].data = hdu[j+1].data[ymin:ymax, xmin:xmax]\n", " hdu[j+1].header['NAXIS1'] = hdu[j+1].data.shape[1]\n", " hdu[j+1].header['NAXIS2'] = hdu[j+1].data.shape[0]\n", " \n", " # Update primary header array size information\n", " hdu[0].header['SIZAXIS1'] = hdu[j+1].data.shape[1]\n", " hdu[0].header['SIZAXIS2'] = hdu[j+1].data.shape[0]\n", " \n", " # Save changes and close\n", " hdu.close()\n", "\n", "# Set drizzle parameters\n", "fs = np.min(np.array(ps)) # Final pixel scale (arcsec) of the output image, native pixel scale: ~0.024742\n", "fp = 1.0 # Fraction by which to shrink the input pixels prior to drizzling onto output grid\n", "\n", "# Drizzle images together\n", "ad.AstroDrizzle(files, static=False, skysub=False, driz_separate=False, median=False, blot=False, \n", " driz_cr=False, driz_combine=True, clean=True, build=False, preserve=False, \n", " final_scale=fs, final_pixfrac=fp, final_wht_type='ERR', output='fuv_drz')" ] }, { "cell_type": "markdown", "id": "bc11f321", "metadata": {}, "source": [ "## About this Notebook\n", "\n", "**Author:** Laura Prichard, Staff Scientist II, STIS Team.\n", "\n", "**Written**: 2022-03-22\n", "\n", "For questions on using the DrizzlePac package with STIS data, contact the [HST Help Desk](https://stsci.service-now.com/hst) (help@stsci.edu).\n", "\n", "## Citations\n", "\n", "- [The DrizzlePac Handbook](https://www.stsci.edu/files/live/sites/www/files/home/scientific-community/software/drizzlepac/_documents/drizzlepac-handbook.pdf): Hoffmann, S. L., Mack, J., et al., 2021, “The DrizzlePac Handbook”, Version, (Baltimore: STScI).\n", "- [STIS Instrument Handbook](https://hst-docs.stsci.edu/stisihb): Prichard, L., Welty, D. and Jones, A., et al. 2022 \"STIS Instrument Handbook,\" Version 21.0, (Baltimore: STScI).\n", "- [STIS Data Handbook](https://hst-docs.stsci.edu/stisdhb): Sohn, S. T., et al., 2019, “STIS Data Handbook”, Version 7.0, (Baltimore: STScI).\n", "- See [citing `astropy`](https://www.astropy.org/acknowledging.html)\n", "_____________\n", "[Top of Page](#top)\n", "\"Space " ] }, { "cell_type": "code", "execution_count": null, "id": "6dc668b4", "metadata": {}, "outputs": [], "source": [] } ], "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.11.8" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": true, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 5 }