diff --git a/notebooks/DrizzlePac/README.md b/notebooks/DrizzlePac/README.md index 99fa42750..5f1d9c8b7 100644 --- a/notebooks/DrizzlePac/README.md +++ b/notebooks/DrizzlePac/README.md @@ -19,6 +19,7 @@ Drizzling Features: - [Creating HST mosaics observed with multiple detectors](https://spacetelescope.github.io/hst_notebooks/notebooks/DrizzlePac/align_mosaics/align_mosaics.html) - [Using Sky Matching features for HST mosaics](https://spacetelescope.github.io/hst_notebooks/notebooks/DrizzlePac/sky_matching/sky_matching.html) - [Masking satellite trails prior to drizzling](https://spacetelescope.github.io/hst_notebooks/notebooks/DrizzlePac/mask_satellite/mask_satellite.html) +- [Drizzling new WFPC2 flt data products](https://spacetelescope.github.io/hst_notebooks/notebooks/DrizzlePac/drizzle_wfpc2/drizzle_wfpc2.ipynb) For more information, see the [DrizzlePac Handbook](https://hst-docs.stsci.edu/drizzpac) and the [readthedocs](https://drizzlepac.readthedocs.io/en/latest/) software documentation. For additional assistance with DrizzlePac tools, users may submit a ticket to the [STScI Help Desk](https://stsci.service-now.com/hst?id=hst_index) and should select the DrizzlePac category. diff --git a/notebooks/DrizzlePac/drizzle_wfpc2/drizzle_wfpc2.ipynb b/notebooks/DrizzlePac/drizzle_wfpc2/drizzle_wfpc2.ipynb index 5666f575c..c66e5d922 100644 --- a/notebooks/DrizzlePac/drizzle_wfpc2/drizzle_wfpc2.ipynb +++ b/notebooks/DrizzlePac/drizzle_wfpc2/drizzle_wfpc2.ipynb @@ -4,117 +4,199 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Drizzling WFPC2 Images to use a Single Zeropoint" + "\n", + "# Drizzling new WFPC2 FLT data products with chip-normalization \"Space\n", + "***\n", + "\n", + "
This notebook requires creating and activating a virtual environment using the requirements file in this notebook's repository. Please also review the README file before using the notebook.
\n", + "\n", + "\n", + "## Table of Contents\n", + "  **[1. Introduction](#intro)
**\n", + "  **[2. Download the Data](#Data)**
\n", + "  **[3. Update the WCS](#WCS)**
\n", + "  **[4. Align the Images](#align)**
\n", + "  **[5. Drizzle-combine the images](#drizzle)**
\n", + "  **[Conclusions](#conclusions)**
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "
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.
" + "## 1. Introduction \n", + "[Table of Contents](#toc)\n", + "\n", + "This notebook shows how to work with a new type of WFPC2 calibrated data product in MAST. The new data products combine the previously separate files containing the science array `c0m.fits`and data quality array `c1m.fits` into a new file with suffix `flt.fits`, similar to calibrated images for ACS and WFC3. \n", + "\n", + "These `flt` files have now been corrected for differences in the inverse sensitivity the science arrays of each chip using the software '[photeq](https://drizzlepac.readthedocs.io/en/deployment/photeq.html)' so that a single PHOTFLAM (or PHOTFNU) value may be used for photometry. \n", + "\n", + "The new `flt` files include four key differences from the older style `c0m`, `c1m` WFPC2 files:\n", + " \n", + " 1.) The SCI & DQ arrays are merged into a single file with multiple extensions.\n", + " 2.) The SCI arrays are corrected for different 'inverse sensitivities' across chips.\n", + " 3.) The SCI arrays from different chips may now be drizzled together in a mosaic.\n", + " 4.) The SCI header WCS includes improved absolute astromentry (i.e. WCSNAME='*GAIA*').\n", + "\n", + " This example uses WFPC2 observations of Omega Centauri (NGC 5139) in the F555W filter taken from CAL proposal [8447](http://www.stsci.edu/cgi-bin/get-proposal-info?id=8447&observatory=HST). The two exposures have been dithered to place the same stars on the WF2 and WF4 chips to measure relative photometry across the two chips. After using `TweakReg` to realign the images, we show how to combine the images with `AstroDrizzle`. We inspect the resulting science and weight images and compare the results with those derived from the older MAST data products.\n", + "\n", + "**Acronyms:**\n", + "- Hubble Space Telescope (HST)\n", + "- Wide Field Camera 3 (WFC3)\n", + "- Advanced Camera for Surveys (ACS)\n", + "- Wide Field and Planetary Camera 2 (WFPC2)\n", + "- Mikulski Archive for Space Telescopes (MAST)\n", + "- Flexible Image Transport System (FITS)\n", + "- World Coordinate System (WCS)\n", + "- Hubble Advanced Products (HAP)\n", + "- Single Visit Mosaic (SVM)\n", + "- Planetary Chip (PC)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - " ## Table of Contents\n", - " [Introduction](#intro)
\n", - " \n", - " [1. Download the Data](#Data)
\n", - " \n", - " [2. Update the WCS](#WCS)
\n", - " - [Backup an Image](#backup)
\n", - " \n", - " [3. Align the Image](#align)
\n", - " [4. Equalize the chip sensitivities](#equalize)
\n", - " [5. Drizzle-combine the images](#drizzle)
\n", - " [6. Illustration of the Effects of Sensitivity Variation on Drizzling](#example)
\n", - " - [Simulate an Image of Constant \"Blank Sky\" Background](#constant)
\n", - " - [Drizzle the Original Simulated Image and the Equalized Image](#drizzling)
\n", - " - [Display The Results of the Simulation Side-by-Side](#display)
\n" + "\n", + "### Import Packages\n", + "[Table of Contents](#toc)\n", + "\n", + "***\n", + "\n", + "The following Python packages are required to run the Jupyter Notebook:\n", + " - [**os**](https://docs.python.org/3/library/os.html) - change and make directories\n", + " - [**glob**](https://docs.python.org/3/library/glob.html) - gather lists of filenames\n", + " - [**shutil**](https://docs.python.org/3/library/shutil.html) - copy files between folders\n", + " - [**numpy**](https://numpy.org) - math and array functions\n", + " - [**matplotlib**](https://matplotlib.org) - create graphics\n", + " - [pyplot](https://matplotlib.org/stable/tutorials/pyplot.html) - make figures and graphics\n", + " - [**astroquery**](https://astroquery.readthedocs.io/en/latest/) - download astronomical data\n", + " - [mast.Observations](https://astroquery.readthedocs.io/en/latest/mast/mast_obsquery.html) - query MAST database\n", + " - [**astropy**](https://www.astropy.org) - model fitting and file handling\n", + " - [io.fits](https://docs.astropy.org/en/stable/io/fits/) - import FITS files\n", + " - [table.Table](https://docs.astropy.org/en/stable/api/astropy.table.QTable.html) - tables without physical units\n", + " - [visualization.LogStretch](https://docs.astropy.org/en/stable/api/astropy.visualization.LogStretch.html#astropy.visualization.LogStretch) - display image normalization\n", + " - [visualization.ImageNormalize](https://docs.astropy.org/en/stable/api/astropy.visualization.ImageNormalize.html#astropy.visualization.ImageNormalize) - display image normalization\n", + " - [**ipython**](https://ipython.org) - cell formatting and interactives\n", + " - [display.Image](https://ipython.readthedocs.io/en/stable/api/generated/IPython.display.html#IPython.display.Image) - display image files\n", + " - [display.clear_output](http://ipython.org/ipython-doc/dev/api/generated/IPython.display.html#IPython.display.clear_output) - clear text from cells\n", + " - [**photutils**](https://photutils.readthedocs.io/en/stable/index.html) - aperture photometry tools\n", + " - [photutils.aperture](https://photutils.readthedocs.io/en/stable/aperture.html) - aperture photometry tools\n", + " - [aperture_photometry](https://photutils.readthedocs.io/en/stable/api/photutils.aperture.aperture_photometry.html#photutils.aperture.aperture_photometry) - measure aperture photometry\n", + " - [CircularAperture](https://photutils.readthedocs.io/en/stable/api/photutils.aperture.CircularAperture.html) - measure photometry in circle\n", + " - [CircularAnnulus](https://photutils.readthedocs.io/en/stable/api/photutils.aperture.CircularAnnulus.html) - measure photometry in annulus\n", + " - [ApertureStats](https://photutils.readthedocs.io/en/stable/api/photutils.aperture.ApertureStats.html) - calculate aperture statistics\n", + " - [**drizzlepac**](https://www.stsci.edu/scientific-community/software/drizzlepac) - combine HST images into mosaics\n", + " - [astrodrizzle](https://drizzlepac.readthedocs.io/en/deployment/astrodrizzle.html) - combine exposures by drizzling\n", + " - [tweakreg](https://drizzlepac.readthedocs.io/en/deployment/tweakreg.html) - align exposures to a common WCS\n", + " - [**stwcs**](https://stwcs.readthedocs.io/en/latest/) - WCS based distortion models and transformations\n", + " - [updatewcs](https://stwcs.readthedocs.io/en/stable/updatewcs.html) - modify WCS in HST file headers" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import glob\n", + "import shutil\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from astroquery.mast import Observations\n", + "from astropy.io import fits\n", + "from astropy.table import Table\n", + "from astropy.visualization import LogStretch, ImageNormalize\n", + "from IPython.display import Image, clear_output\n", + "from photutils.aperture import CircularAnnulus, CircularAperture, ApertureStats, aperture_photometry\n", + "from drizzlepac import tweakreg, astrodrizzle\n", + "from stwcs.updatewcs import updatewcs\n", + "\n", + "%matplotlib inline\n", + "%config InlineBackend.figure_format = 'retina'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Introduction \n", + "## 2. Download the Data \n", + "[Table of Contents](#toc)\n", + " \n", + "The WFPC2 data are downloaded using the `astroquery` API to access the [MAST](http://archive.stsci.edu) archive. The `astroquery.mast` [documentation](http://astroquery.readthedocs.io/en/latest/mast/mast.html) has more examples for how to find and download data from MAST.\n", + "\n", + "MAST queries may be done using `query_criteria`, where we specify:
\n", "\n", - "Extra care must be taken when using `AstroDrizzle` to combine observations from detectors comprised of multiple chips of varying sensitivity. `AstroDrizzle` works with calibrated images in units of counts (electrons or Data Numbers) or count rates and not in units of flux. It assumes that all input frames can be converted to physical flux units using a single inverse-sensitivity factor, recorded in the FITS image headers as `PHOTFLAM`, and the output drizzled product simply copies the `PHOTFLAM` keyword value from the first input image. When this occurs, the inverse-sensitivity will vary across the final drizzled product, and users will need to keep track of which sources fell on which chip when doing photometry. Moreover, varying detector sensitivities will affect the cosmic-ray rejection algorithm used by `AstroDrizzle`, and this may result in the misidentification of some good pixels as cosmic rays. \n", + "    $\\rightarrow$ obs_id, proposal_id, and filters \n", "\n", - "This is a typical situation when drizzle-combining images from HST instruments with different chip sensitivities, e.g. Wide Field and Planetary Camera 2 (WFPC2). For more detail, see the section on [Gain Variation](http://www.stsci.edu/instruments/wfpc2/Wfpc2_dhb/wfpc2_ch53.html) under 'Position-Dependent Photometric Corrections' in the WFPC2 Data Handbook. As a result, each of the four chips requires a [unique PHOTFLAM](http://www.stsci.edu/instruments/wfpc2/Wfpc2_dhb/wfpc2_ch52.html#1933986) header keyword value. A similar situation may occur when drizzle-combining observations taken over a span of several years as detector's sensitivity declines over time, see e.g. [ACS ISR 2016-03](https://doi.org/10.3847/0004-6256/152/3/60).\n", + "MAST data products may be downloaded by using `download_products`, where we specify:
\n", "\n", - "One approach is to rescale the input data so that `AstroDrizzle` can properly assume the images/chips have the same sensitivity; that is, a single `PHOTFLAM` value can be used to convert re-scaled image counts (or count-rates) to physical _integrated_ flux units. The `photeq` task in `Drizzlepac` automates this image intensity rescaling to a single inverse-sensitivity factor `PHOTFLAM`.\n", + "    $\\rightarrow$ data_prod = original `c0m`, `c1m` files and new `flt` calibrated images and their corresponding `drw` drizzled files \n", "\n", - "In this example notebook, archival WFPC2 images are used to demonstrate advanced reprocessing using `TweakReg` and `AstroDrizzle` for alignment and image combination. The notebook is based on a prior WFPC2 [example](https://www.stsci.edu/files/live/sites/www/files/home/scientific-community/software/drizzlepac/examples/_documents/DrizzlePac_EX7.pdf) but includes additional information about equalizing the chip sensitivities prior to combining. \n", + "    $\\rightarrow$ data_type = standard products `CALWFPC2` or Hubble Advanced Products: 'Single Visit Mosaics' `HAP-SVM`\n", "\n", - "**NOTE:** It is important to note that `photeq` only adjusts image counts so that _integrated_ physical fluxes can be obtained using a single `PHOTFLAM`. It does nothing to account for different throughtputs at different wavelengths." + "
Depending on your connection speed this cell may take a few minutes to execute.
" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ - "import shutil\n", - "import glob\n", - "import os\n", + "obs_ids = ['U5JX0108R', 'U5JX010HR']\n", + "filts = ['F555W']\n", "\n", - "import matplotlib.pyplot as plt\n", - "from astropy.io import fits\n", + "obsTable = Observations.query_criteria(obs_id=obs_ids, filters=filts)\n", + "products = Observations.get_product_list(obsTable)\n", "\n", - "from astroquery.mast import Observations\n", - "from stwcs.updatewcs import updatewcs\n", - "from drizzlepac import tweakreg, astrodrizzle, photeq\n", + "data_prod = ['C0M', 'C1M', 'FLT', 'DRW'] \n", + "data_type = ['CALWFPC2'] # Options: ['CALWFPC2', 'HAP-SVM']\n", + " \n", + "Observations.download_products(products, productSubGroupDescription=data_prod, project=data_type)\n", "\n", - "# ONLY needed for the simulation section:\n", - "import numpy as np\n", - "from stwcs.wcsutil import HSTWCS\n", - "from drizzlepac.wfpc2Data import WFPC2_GAINS\n", - "\n", - "%matplotlib inline" + "clear_output()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## 1. Download the Data \n", + "The MAST pipeline uses the new `flt` files to create two new types of WFPC2 drizzled products. Unlike the original `u*_drz.fits` files, which were drizzled to an output scale of the WF chips (0.0996\"/pix), the new drizzled files `u*_drw.fits` are drizzled to an output scale of the PC chip (0.0455\"/pix), where the new `drw` suffix reflects the improved WCS. \n", "\n", - "This example uses WFPC2 observations of Messier 2 in the F814W filter. The data come from GO proposal [11100](http://www.stsci.edu/cgi-bin/get-proposal-info?id=11100&observatory=HST) _\"Two new 'bullets' for MOND: revealing the properties of dark matter in massive merging clusters\"_. Four images were acquired using a 4-pt dither box pattern, followed by two images offset with a dither-line pattern. \n", + "A second drizzled product `hst_*drz.fits` combines all images in the same filter and may include improved relative alignment across filters in the same visit. These are referred to as Hubble Advanced Products and are listed below, but are beyond the scope of this notebook. For more details, see the [HAP Webpage](https://outerspace.stsci.edu/pages/viewpage.action?spaceKey=HAdP&title=Improvements+in+HST+Astrometry) and the [Drizzlepac Handbook](https://hst-docs.stsci.edu/drizzpac).\n", + " \n", "\n", - "The data are downloaded using the `astroquery` API to access the [MAST](http://archive.stsci.edu) archive. The `astroquery.mast` [documentation](http://astroquery.readthedocs.io/en/latest/mast/mast.html) has more examples for how to find and download data from MAST." + " Drizzled Image Dimensions Scale Input Updated\n", + " (\"/pix) File WCS?\n", + " ______________________________________________________________________________________\n", + " u5jx010hr_drz.fits [1515, 1495] 0.0996 c?m.fits N\n", + " u5jx010hr_drw.fits [4171, 4143] 0.0455 flt.fits Y\n", + " hst_8447_01_wfpc2_pc_f555w_u5jx01_drz.fits [6258, 6215] 0.0455 flt.fits Y " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Move files to the current working directory" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ - "# Retrieve the observation information.\n", - "if os.path.isdir('mastDownload'):\n", - " shutil.rmtree('mastDownload')\n", - "obs_table = Observations.query_criteria(obs_id='ua0605*', filters='F814W', obstype='ALL')\n", - "products = Observations.get_product_list(obs_table)\n", - "\n", - "# Download only the ua0605*_c0m.fits and ua0605*_c1m.fits (DQ) images:\n", - "Observations.download_products(products, mrp_only=False, productSubGroupDescription=['C0M', 'C1M'], extension='fits')\n", - "\n", - "# Move the files from the mastDownload directory to the current working\n", - "# directory and make a backup of the files.\n", - "fits_files = glob.glob('mastDownload/HST/ua*/ua*c?m.fits')\n", - "for fil in fits_files:\n", - " base_name = os.path.basename(fil)\n", + "fits_files = glob.glob('mastDownload/HST/*/*.fits')\n", + "for f in fits_files:\n", + " base_name = os.path.basename(f)\n", " if os.path.isfile(base_name):\n", " os.remove(base_name)\n", - " shutil.move(fil, '.')\n", - " \n", - "# Delete the mastDownload directory and all subdirectories it contains.\n", + " shutil.move(f, '.') \n", "shutil.rmtree('mastDownload')" ] }, @@ -122,54 +204,88 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## 2. Update the WCS \n", + "### Inspect the File Structure\n", + "\n", + "While the `flt.fits` files use the same notation for different chips, `EXT=('SCI',1)` the order of the FITS extensions (`'SCI', 'ERR', 'DQ'`) is different for WFC3/UVIS and WFPC2. For example:\n", "\n", - "WFPC2 images downloaded from the archive contain World Coordinate System (WCS) information based on an older-style description of image distortions. Before these images can be processed with `drizzlepac`, their WCS must be converted to a new format. This can be achieved using `updatewcs()` function from the `stwcs` package. More details may be found in the DrizzlePac Handbook v1.0, chapter 4.4.5 under section:['Making WFPC2 Images Compatible with AstroDrizzle'](https://www.stsci.edu/files/live/sites/www/files/home/scientific-community/software/drizzlepac/_documents/drizzlepac-handbook-v1.pdf). Note that `updatewcs` is no longer a parameter in AstroDrizzle or TweakReg and must be run separately before processing the data. \n", + " uvis_flt.fits[1] = uvis_flt.fits['SCI',1] --> UVIS2, SCI\n", + " uvis_flt.fits[2] = uvis_flt.fits['ERR',1] --> UVIS2, ERR\n", + " uvis_flt.fits[3] = uvis_flt.fits['DQ',1] --> UVIS2, DQ\n", + " uvis_flt.fits[4] = uvis_flt.fits['SCI',2] --> UVIS1, SCI\n", + " uvis_flt.fits[5] = uvis_flt.fits['ERR',2] --> UVIS1, ERR\n", + " uvis_flt.fits[6] = uvis_flt.fits['DQ',2] --> UVIS1, DQ\n", "\n", - "First we download the reference files from the CRDS website. See the initialization notebook in this repository for more information. " + " wfpc2_flt.fits[1] = wfpc2_flt.fits['SCI',1) --> PC1, SCI\n", + " wfpc2_flt.fits[2] = wfpc2_flt.fits['SCI',2] --> WF2, SCI\n", + " wfpc2_flt.fits[3] = wfpc2_flt.fits['SCI',3] --> WF3, SCI\n", + " wfpc2_flt.fits[4] = wfpc2_flt.fits['SCI',4] --> WF4, SCI\n", + " wfpc2_flt.fits[5] = wfpc2_flt.fits['DQ',1] --> PC1, DQ\n", + " wfpc2_flt.fits[6] = wfpc2_flt.fits['ERR',1] --> PC1, ERR \n", + " wfpc2_flt.fits[7] = wfpc2_flt.fits['DQ',2] --> WF2, DQ\n", + " wfpc2_flt.fits[8] = wfpc2_flt.fits['ERR',2] --> WF2, ERR\n", + " wfpc2_flt.fits[9] = wfpc2_flt.fits['DQ',3] --> WF3, DQ\n", + " wfpc2_flt.fits[10] = wfpc2_flt.fits['ERR',3] --> WF3, ERR\n", + " wfpc2_flt.fits[11] = wfpc2_flt.fits['DQ',4] --> WF4, DQ\n", + " wfpc2_flt.fits[12] = wfpc2_flt.fits['ERR',4] --> WF4, ERR" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ - "os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'\n", - "os.environ['CRDS_PATH'] = os.path.abspath(os.path.join('.', 'reference_files'))\n", - "os.system('crds bestrefs --files ua0605*_c0m.fits --sync-references=1 --update-bestrefs')\n", - "os.environ['uref'] = os.path.abspath(os.path.join('.', 'reference_files', 'references', 'hst', 'wfpc2')) + os.path.sep" + "fits.info('u5jx010hr_c0m.fits')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fits.info('u5jx010hr_flt.fits')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "**NOTE:** This next step may raise warnings because the Astrometry database is in progress and currently does not cover the entire sky. Please ignore these warnings. The WCS will still be updated. " + "The drizzled products, on the other hand, contain 3 data extensions: `SCI`, `WHT`, and `CTX`. " ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ - "updatewcs('ua*c0m.fits', use_db=True)" + "fits.info('u5jx010hr_drw.fits')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Backup an Image \n", + "### Compare photometric keywords from the handbook with the image headers\n", + " \n", + "`PHOTFLAM` 'inverse sensitivity' values are listed in Table 5.1 of the WFPC2 DATA HANDBOOK Section 5.1 [Photometric Zeropoint](http://www.stsci.edu/instruments/wfpc2/Wfpc2_dhb/wfpc2_ch53.html) for the gain 7 setting (which is the default setting for most science programs). The `PHOTFLAM` values for gain 14 can be obtained by multiplying by the gain ratios: 1.987 (PC1), 2.003 (WF2), 2.006 (WF3), and 1.955 (WF4), where the values are from [Holtzman et al. (1995)](https://ui.adsabs.harvard.edu/abs/1995PASP..107.1065H/abstract). \n", "\n", - "In a later section we will generate simulated data to illustrate the effects of drizzling WFPC2 images without sensitivity equalization. For that purpose we will need a copy of an original image that has the original inverse-sensitivity values (`PHOTFLAM`) in their headers. Here we create a backup copy of the first image.\n", + " Filter PC1 WF2 WF3 WF4\n", + " ----- --------- --------- --------- ---------\n", + " f555w \t 3.483e-18 3.396e-18 3.439e-18 3.507e-18 \n", "\n", - "**NOTE:** This step is needed for illustration purpose in this notebook only. It is not needed when processing data." + "The `PHOTFLAM` values include differences in the gain (e- DN-1) between chips, i.e. Table 4.2 in Section 4.4: \"Read Noise and Gain Settings\" in the [2008 WFPC2 Instrument Handbook](https://www.stsci.edu/files/live/sites/www/files/home/hst/instrumentation/legacy/wfpc2/_documents/wfpc2_ihb.pdf), reproduced below:\n", + "\n", + " Gain PC1 WF2 WF3 WF4\n", + " ----- --------- --------- --------- ---------\n", + " \"7\" 7.12 ± 0.41 7.12 ± 0.41 6.90 ± 0.32 7.10 ± 0.39\n", + " \"15\" 13.99 ± 0.63 14.50 ± 0.77 13.95 ± 0.63 13.95 ± 0.70\n", + " \n", + "\n", + "While WFC3 and ACS flt.fits data products are in units of electrons (WFC3/UVIS, ACS/WFC) or electrons/sec (WFC3/IR), WFPC2 images are in units of 'Data Number' (DN) or 'COUNTS' as reflected in the BUNIT keyword.
\n", + " \n", + "Here we inspect a subset of photometry keywords in the `c0m` and `flt` files:" ] }, { @@ -178,45 +294,167 @@ "metadata": {}, "outputs": [], "source": [ - "orig_image = glob.glob('ua*c0m.fits')[0]\n", - "backup_image = 'simulation.fits'\n", - "if os.path.isfile(backup_image):\n", - " os.remove(backup_image)\n", - "shutil.copy2(orig_image, backup_image)" + "with fits.open('u5jx0108r_c0m.fits') as h:\n", + " phot1 = h[1].header['PHOTFLAM']\n", + " phot2 = h[2].header['PHOTFLAM']\n", + " phot3 = h[3].header['PHOTFLAM']\n", + " phot4 = h[4].header['PHOTFLAM']\n", + " gain = h[0].header['ATODGAIN']\n", + " bunit = h[1].header['BUNIT']\n", + "\n", + "print('PHOTFLAM for PC1: ', phot1, ' Factor Relative to WF3: ', round(phot1/phot3, 4))\n", + "print('PHOTFLAM for WF2: ', phot2, ' Factor Relative to WF3: ', round(phot2/phot3, 4))\n", + "print('PHOTFLAM for WF3: ', phot3, ' Factor Relative to WF3: ', round(phot3/phot3, 4))\n", + "print('PHOTFLAM for WF4: ', phot4, ' Factor Relative to WF3: ', round(phot4/phot3, 4))\n", + "print('Gain Value is: ', gain) \n", + "print('Units are: ', bunit) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## 3. Align the Images \n", + "Note the unique PHOTFLAM values for each chip in the original `c0m` data products (above). \n", "\n", - "Due to small pointing errors, the image header WCS typically needs to be updated in order to achieve the best drizzle-combined products. The expected pointing accuracy for various observing scenerios is summarized in the DrizzlePac Handbook [Chapter 4.4](https://hst-docs.stsci.edu/drizzpac/chapter-4-astrometric-information-in-the-header/4-4-hst-pointing-accuracy-and-stability). Input images must first be aligned so that when the coordinates of a given object (in detector space) are converted to sky coordinates (using the WCS), that object's sky coordinates must be approximately equal in each frame. \n", + "In the new `flt` data products, the PHOTFLAM values are all set to the value for WF3 (below)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "with fits.open('u5jx0108r_flt.fits') as h:\n", + " phot1 = h[1].header['PHOTFLAM']\n", + " phot2 = h[2].header['PHOTFLAM']\n", + " phot3 = h[3].header['PHOTFLAM']\n", + " phot4 = h[4].header['PHOTFLAM']\n", + " gain = h[0].header['ATODGAIN']\n", + " bunit = h[1].header['BUNIT']\n", "\n", - "The `DrizzlePac` task `TweakReg` may be used to correct for any errors in the image header WCS. First, `TweakReg` finds sources in each image, matches sources in common across images, and finds a separate linear transformation to align each image. `TweakReg` then computes a new WCS for each image based on this linear transformation.\n", + "print('PHOTFLAM for PC1: ', phot1, ' Factor Relative to WF3: ', round(phot1/phot3, 4))\n", + "print('PHOTFLAM for WF2: ', phot2, ' Factor Relative to WF3: ', round(phot2/phot3, 4))\n", + "print('PHOTFLAM for WF3: ', phot3, ' Factor Relative to WF3: ', round(phot3/phot3, 4))\n", + "print('PHOTFLAM for WF4: ', phot4, ' Factor Relative to WF3: ', round(phot4/phot3, 4))\n", + "print('Gain Value is: ', gain) \n", + "print('Units are: ', bunit) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The new `drw` data products carry the WF3 PHOTFLAM values. While the BUNIT keyword says 'counts', it should say 'counts/sec'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "with fits.open('u5jx010hr_drw.fits') as h:\n", + " phot1 = h[1].header['PHOTFLAM']\n", + " gain = h[0].header['ATODGAIN']\n", + " bunit = h[1].header['BUNIT']\n", + "\n", + "print('PHOTFLAM for Combined image', phot1, ' Factor Relative to WF3: ', round(phot1/phot3, 4))\n", + "print('Gain Value is: ', gain) \n", + "print('Units are: ', bunit) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Because WF3 is the chip with the most accurate flux calibration, this was selected as the 'reference' chip and the inverse sensitivities in the new `flt` files have been equalized to that value. For example:\n", + "\n", + " photeq.photeq(files='*_flt.fits', ref_phot_ext=3) # or ('sci',3)\n", + " \n", + "\n", + "Note that `photeq` adjusts the pixels values in the SCI array so that photometry can be obtained using a single `PHOTFLAM`. The ratio of the count rates in the chips will therefore be equal to the ratio of the PHOTFLAM values. \n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Compare the equation for computing STMAG using the new and the old workflow: \n", "\n", - "Here we show a basic image alignment procedure. For a more detailed illustration of image alignment, please refer to other example notebooks such as the mosaic example in this repository." + "NEW: Use `flt` files \n", + " \n", + "Drizzle the 'flt.fits' files and do photometry for all chips at once; \n", + "\n", + " STMAG = -21.1 -2.5*log (countrate * PHOTFLAM) \n", + "\n", + "OLD: Use `c0m` files \n", + " \n", + " option A.) Drizzle the 'c0m.fits' files and do photometry for each chip separately; \n", + " \n", + " STMAG = -21.1 -2.5*log (countrate * PHOTFLAM_1,2,3,4) \n", + " \n", + " option B.) Run [phot_eq](https://drizzlepac.readthedocs.io/en/deployment/photeq.html) on the 'c0m.fits' files prior to drizzling, normalizing to the WF3 chip. Use the first equation with a single PHOTFLAM value which equal to PHOTFLAM for WF3 in the second equation.\n", + " \n", + "With the older `c0m` data products, caution must be taken when using `AstroDrizzle` to combine chips of different sensitivity. `AstroDrizzle` works with calibrated images in units of counts (electrons or Data Numbers) or count rates and not in units of flux. It assumes that all input frames can be converted to physical flux units using a single inverse-sensitivity value, and the output drizzled product simply copies the `PHOTFLAM` keyword value from the first input image. When this occurs, the inverse-sensitivity will vary across the final drizzled product, and users will need to keep track of which sources fell on which chip when doing photometry. Moreover, varying detector sensitivities can affect the cosmic-ray rejection algorithm used by `AstroDrizzle`, and this may result in the misidentification of some good pixels as cosmic rays. Using the new FLT style files solves these issues by normalizing the chips to a single `PHOTFLAM` value." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Check image header data \n", + "\n", + "In this section, we inspect keywords in the image headers which are useful for understanding the data. Note the `PHOTMODE` keyword in column 6 below, which lists the CCD gain. Additionally, we inspect the keywords `WCSNAME` and `NMATCHES` to determine whether the images have improved astrometry." ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ - "tweakreg.TweakReg('ua*c0m.fits', updatehdr=True, reusename=True, interactive=False,\n", - " conv_width=3.0, threshold=300.0, peakmin=100, peakmax=10000)" + "paths = sorted(glob.glob('*flt.fits'))\n", + "data = []\n", + "keywords_ext0 = [\"ROOTNAME\", \"TARGNAME\", \"FILTNAM1\", \"EXPTIME\"]\n", + "keywords_ext1 = [\"ORIENTAT\", \"PHOTMODE\", \"PHOTFLAM\", \"WCSNAME\", \"NMATCHES\"]\n", + "\n", + "for path in paths:\n", + " path_data = []\n", + " for keyword in keywords_ext0:\n", + " path_data.append(fits.getval(path, keyword, ext=0))\n", + " for keyword in keywords_ext1:\n", + " path_data.append(fits.getval(path, 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', 'f8', 'f8', 'str', 'f8', 'str', 'i8'])\n", + "table['EXPTIME'].format = '7.1f' \n", + "table['ORIENTAT'].format = '7.2f'\n", + "table['PHOTFLAM'].format = '7.6g' \n", + "table" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## 4. Equalize the chip sensitivities \n", + "Here we see that the images have been aligned to Gaia eDR3 with ~ 100 matched objects in the HST data." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3. Update the WCS \n", + "[Table of Contents](#toc)\n", + "\n", + "Unlike the new `flt` products, the older WFPC2 `c0m` files contain World Coordinate System (WCS) information based on an older-style description of image distortions. Before these images can be processed with `drizzlepac`, their WCS must be converted to a new format. This can be achieved using `updatewcs()` function from the `stwcs` package. \n", + "\n", + "More details may be found in the DrizzlePac Handbook, chapter 4.4.5 under section: ['Making WFPC2 Images Compatible with AstroDrizzle'](https://www.stsci.edu/files/live/sites/www/files/home/scientific-community/software/drizzlepac/_documents/drizzlepac-handbook-v1.pdf). NOTE: Running `updatewcs` is not necessary when working with `flt` files. \n", "\n", - "This step adjusts image data values so that all images and chips appear (to `AstroDrizzle`) to have a single inverse sensitivity (`PHOTFLAM`). This can be achieved using the `photeq` task in `Drizzlepac`. This task adjusts image data so that when these data are multiplied by the same single `PHOTFLAM` value, the correct flux is obtained." + "First we download the necessary reference files from the CRDS website." ] }, { @@ -227,30 +465,62 @@ }, "outputs": [], "source": [ - "photeq.photeq(files='ua*_c0m.fits', ref_phot_ext=3, readonly=False)" + "os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'\n", + "os.environ['CRDS_PATH'] = os.path.abspath(os.path.join('.', 'reference_files'))\n", + "os.system('crds bestrefs --files u5jx01*_c0m.fits --sync-references=1 --update-bestrefs')\n", + "os.environ['uref'] = os.path.abspath(os.path.join('.', 'reference_files', 'references', 'hst', 'wfpc2')) + os.path.sep\n", + "\n", + "clear_output()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "In the above command, we instruct `photeq` to \"equalize\" all chips of all input images using the `PHOTFLAM` for the `WF3` chip (`ref_phot_ext=3`), using the first image as a reference. This reference `PHOTFLAM` value is reported in the log file:\n", - "\n", - "```\n", - "REFERENCE VALUE FROM FILE: 'ua060502m_c0m.fits['SCI',1]'\n", - "REFERENCE 'PHOTFLAM' VALUE IS: 2.507987e-18\n", - "```\n", - "\n", - "Upon the completion, `photeq` will not only adjust image data but also update `PHOTFLAM` values for all chips to this specific reference value." + "**NOTE:** This next step may raise warnings which may be ignored. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "updatewcs('u*c0m.fits', use_db=False)\n", + "clear_output()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## 5. Drizzle-combine the images \n", + "While the new `flt` images have already been aligned to Gaia, we want to directly compare (ratio) the two sets of drizzled products. Thus we reset the WCS in the `flt` files to use the 'distortion-only' solution. When there are a large number of GAIA matches ($\\geq$ 20), users will want to use the GAIA-based WCS solutions for optimal alignment, especially when there is a large dither across chips, which can introduce a significant rotational offset between exposures (as seen in the `TweakReg` shift files below). In this notebook we focus on photometric differences between the WF2 and WF4 chips, so simply align the exposures using bright stars in common between the two `flt` images." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "updatewcs('u*flt.fits', use_db=False)\n", + "clear_output()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 4. Align the Images \n", + "[Table of Contents](#toc)\n", "\n", - "All four chips are now drizzled together with an output pixel scale set to that of the WF chips:" + "After resetting the WCS, we check for any residual pointing errors between the two dithered exposures. Precise aligment is needed to achieve the best drizzle-combined products. The expected pointing accuracy for various observing scenerios is summarized in the DrizzlePac Handbook [Chapter 4.4](https://hst-docs.stsci.edu/drizzpac/chapter-4-astrometric-information-in-the-header/4-4-hst-pointing-accuracy-and-stability). Input images must first be aligned so that when the coordinates of a given object (in detector space) are converted to sky coordinates (using the WCS), that object's sky coordinates must be approximately equal in each frame. \n", + "\n", + "The `DrizzlePac` task `TweakReg` may be used to correct for any errors in the image header WCS. First, `TweakReg` finds sources in each image, matches sources in common across images, and finds a separate linear transformation to align each image. `TweakReg` then computes a new WCS for each image based on this linear transformation.\n", + "\n", + "Here we show a basic image alignment procedure. For a more detailed illustration of image alignment, please refer to other example notebooks in this repository. Note we use a farily large `searchrad` of 1 arcsec, because we have removed the Gaia WCS and are using the original 'distortion only' solution which can have large errors in the case of large dithers. \n", + "\n", + "First we align the `flt` files and inspect the results before aligning the `c0m` files and verifying that output `shiftfile` values are identical. While testing the alignment, we recommend setting `updatehdr` = `False` to confirm that the alignment results look accurate and then rerunning with the keyword set to `True`." ] }, { @@ -261,20 +531,29 @@ }, "outputs": [], "source": [ + "input_images = sorted(glob.glob('u*flt.fits'))\n", "\n", - "astrodrizzle.AstroDrizzle('ua*c0m.fits',\n", - " preserve=False,\n", - " driz_sep_bits='8,1024',\n", - " driz_sep_wcs=True,\n", - " driz_sep_scale=0.0996,\n", - " combine_type='median',\n", - " driz_cr_snr='5.5 3.5',\n", - " driz_cr_scale='2.0 1.5',\n", - " final_fillval=None,\n", - " final_bits='8,1024',\n", - " final_wcs=True,\n", - " final_scale=0.0996,\n", - " skysub=False)" + "tweakreg.TweakReg(input_images, \n", + " updatehdr=True, \n", + " clean=True,\n", + " reusename=True, \n", + " interactive=False,\n", + " conv_width=3.0, \n", + " threshold=200.0,\n", + " ylimit=1,\n", + " shiftfile=True, \n", + " outshifts='shift_flt.txt',\n", + " searchrad=1,\n", + " tolerance=3,\n", + " minobj=7)\n", + "clear_output()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### If the alignment is unsuccessful, stop the notebook" ] }, { @@ -283,30 +562,48 @@ "metadata": {}, "outputs": [], "source": [ - "drw_data = fits.getdata('final_drw_sci.fits')\n", - "plt.figure(figsize=(10, 10))\n", - "plt.imshow(drw_data, cmap='gray', vmin=-0.1, vmax=0.5, origin='lower')\n", - "plt.show()" + "with open('shift_flt.txt', 'r') as shift:\n", + " for line_number, line in enumerate(shift, start=1):\n", + " if \"nan\" in line:\n", + " raise ValueError('nan found in line {} in shift file'.format(line_number))\n", + " else:\n", + " continue" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## 6. Illustration of the Effects of Sensitivity Variation on Drizzling \n", + "### Inspect the shift file and fit residuals\n", "\n", - "The effect of drizzling images with different detector sensitivies, while tangible, is sometimes difficult to _visualize_ in noisy data, especially when drizzling multiple dithered images that can blur the borders between chips.\n", + "We can look at the shift file to see how well the fit did (or we could open the output png images for more information). The columns are:\n", + "- Filename, X Shift [pixels], Y Shift [pixels], Rotation [degrees], Scale, X RMS [pixels], Y RMS [pixels]\n", "\n", - "In this section we produce a simple simulation of observing a constant intensity \"blank sky\". We then make a copy of this image and apply sensitivity equalization to it. Finally we drizzle both images and compare them side-by-side." + "**Note that there is a small shift in x=1.8 pix and y=-0.4 pix and large rotation of 0.014 degrees between the two images, likely due to the large chip-sized dither.** These residuals are larger than the RMS scatter of 0.15 pixels and are therefore likely to be real. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "shift_flt_table = Table.read('shift_flt.txt', format='ascii.no_header', \n", + " names=['file', 'dx', 'dy', 'rot', 'scale', 'xrms', 'yrms'])\n", + "\n", + "formats = ['.2f', '.2f', '.3f', '.5f', '.2f', '.2f']\n", + "for i, col in enumerate(shift_flt_table.colnames[1:]):\n", + " shift_flt_table[col].format = formats[i]\n", + "shift_flt_table" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Simulate an Image of Constant \"Blank Sky\" Background \n", + "### Inspect the Astrometric residual scatter plots\n", "\n", - "In this simple simulation we assume only Poisson noise." + "Good practice includes inspecting not only the shift file residuals, but also the accompanying diagnostic plots. The residuals in the plot below are clustered around 0.0 with no obvious systematics." ] }, { @@ -315,63 +612,14 @@ "metadata": {}, "outputs": [], "source": [ - "with fits.open(backup_image, mode='update') as h:\n", - " # get chip inverse-sensitivity:\n", - " phot1 = h[1].header['PHOTFLAM']\n", - " phot2 = h[2].header['PHOTFLAM']\n", - " phot3 = h[3].header['PHOTFLAM']\n", - " phot4 = h[4].header['PHOTFLAM']\n", - " \n", - " # get chip WCS:\n", - " w1 = HSTWCS(h, ext=1)\n", - " w2 = HSTWCS(h, ext=2)\n", - " w3 = HSTWCS(h, ext=3)\n", - " w4 = HSTWCS(h, ext=4)\n", - " ref_pscale = w4.idcscale\n", - " \n", - " # get chip gain:\n", - " cmdgain = h[0].header['ATODGAIN']\n", - " gain1 = WFPC2_GAINS[1][cmdgain][0]\n", - " gain2 = WFPC2_GAINS[2][cmdgain][0]\n", - " gain3 = WFPC2_GAINS[3][cmdgain][0]\n", - " gain4 = WFPC2_GAINS[4][cmdgain][0]\n", - " \n", - " # final drizzle scale:\n", - " scale = 0.0996\n", - " \n", - " # simulated sky level (\"true\" sky is constant):\n", - " sky = 10 * phot3\n", - "\n", - " # simulate observed counts assuming only Poisson noise:\n", - " conv1a = gain1 * (w1.idcscale / ref_pscale)**2 / phot1\n", - " conv1b = (gain4 / gain1**2) * (ref_pscale / scale)**2\n", - " conv2a = gain2 * (w2.idcscale / ref_pscale)**2 / phot2\n", - " conv2b = (gain4 / gain2**2) * (ref_pscale / scale)**2\n", - " conv3a = gain3 * (w3.idcscale / ref_pscale)**2 / phot3\n", - " conv3b = (gain4 / gain3**2) * (ref_pscale / scale)**2\n", - " conv4a = gain4 * (w4.idcscale / ref_pscale)**2 / phot4\n", - " conv4b = (1.0 / gain4) * (ref_pscale / scale)**2\n", - "\n", - " h[1].data[:, :] = np.random.poisson(conv1a * sky, h[1].data.shape) * conv1b\n", - " h[2].data[:, :] = np.random.poisson(conv2a * sky, h[2].data.shape) * conv2b\n", - " h[3].data[:, :] = np.random.poisson(conv3a * sky, h[3].data.shape) * conv3b\n", - " h[4].data[:, :] = np.random.poisson(conv4a * sky, h[4].data.shape) * conv4b\n", - "\n", - "# make a copy of this file:\n", - "photeq_image = 'simulation_eq.fits'\n", - "if os.path.isfile(photeq_image):\n", - " os.remove(photeq_image)\n", - "shutil.copy2(backup_image, photeq_image)\n", - "\n", - "# apply equalization to the image copy:\n", - "photeq.photeq(files=photeq_image, ref_phot_ext=3, readonly=False)" + "Image(filename='residuals_u5jx010hr_flt.png', width=500, height=300)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We need to make a simulated DQ array. We used the first science image in the downloaded dataset for the simulation, so we will use its DQ image. " + "Here we align the `c0m` files with the same parameter settings:" ] }, { @@ -382,34 +630,97 @@ }, "outputs": [], "source": [ - "# make an DQ array for every \n", - "hdu = fits.open('ua06050cm_c1m.fits')\n", - "shape = np.shape(hdu[1].data)\n", - "zero_DQ = np.zeros(shape)\n", + "input_images = sorted(glob.glob('u*c0m.fits'))\n", "\n", - "fits_file = fits.HDUList()\n", - "# copy all of its extensions over, including the header\n", - "hdu0 = fits.PrimaryHDU(header=hdu[0].header)\n", - "hdu1 = fits.ImageHDU(hdu[1].data, header=hdu[1].header)\n", - "hdu2 = fits.ImageHDU(hdu[2].data, header=hdu[2].header)\n", - "hdu3 = fits.ImageHDU(hdu[3].data, header=hdu[3].header)\n", - "hdu4 = fits.ImageHDU(hdu[4].data, header=hdu[4].header)\n", + "tweakreg.TweakReg(input_images, \n", + " updatehdr=True, \n", + " clean=True,\n", + " reusename=True, \n", + " interactive=False,\n", + " conv_width=3.0, \n", + " threshold=200.0,\n", + " ylimit=1,\n", + " shiftfile=True, \n", + " outshifts='shift_c0m.txt',\n", + " searchrad=1,\n", + " tolerance=3,\n", + " minobj=7)\n", "\n", - "fits_file.append(hdu0)\n", - "fits_file.append(hdu1)\n", - "fits_file.append(hdu2)\n", - "fits_file.append(hdu3)\n", - "fits_file.append(hdu4)\n", + "clear_output()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, we compare the shift file computed from the `flt` and `c0m` files to confirm that the results are identical." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "shift_c0m_table = Table.read('shift_c0m.txt', format='ascii.no_header', \n", + " names=['file', 'dx', 'dy', 'rot', 'scale', 'xrms', 'yrms'])\n", "\n", - "fits_file.writeto('simula_c1n.fits', overwrite=True)\n", - "fits_file.writeto('simulatio_c1q.fits', overwrite=True)" + "formats = ['.2f', '.2f', '.3f', '.5f', '.2f', '.2f']\n", + "for i, col in enumerate(shift_c0m_table.colnames[1:]):\n", + " shift_c0m_table[col].format = formats[i]\n", + " \n", + "shift_c0m_table" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "shift_flt_table" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Drizzle the Original Simulated Image and the Equalized Image " + "## 5. Drizzle-combine the images \n", + "[Table of Contents](#toc)\n", + "\n", + "All four chips are now drizzled together. For speed in the notebook, we set the output pixel scale to that of the WF chips, but this can be reset to the PC scale if desired, to match the scale of the new `drw` drizzled products. \n", + "\n", + "While we have turned sky subtraction off in the notebook, recommend using `skymethod` = 'localmin' and not 'match' for WFPC2 to avoid issues with overlapping vignetted regions along chip boundaries. See Figure 3.1.2 in the [WFPC2 Instrument Handbook](https://www.stsci.edu/files/live/sites/www/files/home/hst/instrumentation/legacy/wfpc2/_documents/wfpc2_ihb.pdf). For more information drizzling HST data, see the [Drizzlepac Handbook](https://hst-docs.stsci.edu/drizzpac) and the list of commands in the [readthedocs](https://drizzlepac.readthedocs.io/en/latest/drizzlepac_api/astrodrizzle.html).\n", + "\n", + "Users may wish turn on the sky subtraction and cosmic-ray rejection steps by setting `skysub`, `median`, `blot`, and `driz_cr` parameters to `True`. In that case, the recommended parameter settings for WFPC2 are `driz_cr_snr='5.5 3.5'`, and `driz_cr_scale='2.0 1.5'`. For more discussion on drizzling parameters for WFPC2 when using the output scale of the WF chips, see the [Prior Drizzling Example](https://www.stsci.edu/files/live/sites/www/files/home/scientific-community/software/drizzlepac/examples/_documents/DrizzlePac_EX7.pdf)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "astrodrizzle.AstroDrizzle('u*flt.fits', \n", + " output='wfpc2_flt', \n", + " preserve=False, \n", + " build=False, \n", + " context=False, \n", + " skysub=False, \n", + " driz_sep_wcs=True, \n", + " driz_sep_scale=0.0996, \n", + " driz_sep_bits='8,1024', \n", + " driz_sep_fillval=-1, \n", + " median=False, \n", + " blot=False,\n", + " driz_cr=False, \n", + " final_fillval=None, \n", + " final_bits='8,1024', \n", + " final_wcs=True, \n", + " final_scale=0.0996)\n", + "clear_output()" ] }, { @@ -420,50 +731,84 @@ }, "outputs": [], "source": [ - "astrodrizzle.AstroDrizzle(\n", - " backup_image,\n", - " output='nonequalized.fits',\n", - " stepsize=1,\n", - " preserve=False,\n", - " restore=False,\n", - " in_memory=True,\n", - " context=False,\n", - " build=False,\n", - " static=False,\n", - " skysub=False,\n", - " median=False,\n", - " blot=False,\n", - " driz_cr=False,\n", - " final_fillval=None,\n", - " final_bits='',\n", - " final_wcs=True,\n", - " final_scale=scale)\n", + "astrodrizzle.AstroDrizzle('u*c0m.fits', \n", + " output='wfpc2_c0m', \n", + " preserve=False, \n", + " build=False, \n", + " context=False, \n", + " static=True, \n", + " skysub=False, \n", + " driz_sep_wcs=True, \n", + " driz_sep_refimage='u5jx010hr_single_sci.fits', \n", + " driz_sep_bits='8,1024', \n", + " driz_sep_fillval=-1, \n", + " median=False, \n", + " blot=False, \n", + " driz_cr=False, \n", + " final_fillval=None, \n", + " final_bits='8,1024', \n", + " final_wcs=True, \n", + " final_scale=0.0996, \n", + " final_refimage='wfpc2_flt_drw_sci.fits')\n", + "clear_output()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Clean up unnecessary mask files." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fits_files = glob.glob('*ask.fits')\n", + "for f in fits_files:\n", + " base_name = os.path.basename(f)\n", + " if os.path.isfile(base_name):\n", + " os.remove(base_name)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Inspect the single drizzled output frames derived from the dithered `flt` inputs\n", "\n", - "astrodrizzle.AstroDrizzle(\n", - " photeq_image,\n", - " output='equalized.fits',\n", - " stepsize=1,\n", - " preserve=False,\n", - " restore=False,\n", - " in_memory=True,\n", - " context=False,\n", - " build=False,\n", - " static=False,\n", - " skysub=False,\n", - " median=False,\n", - " blot=False,\n", - " driz_cr=False,\n", - " final_fillval=None,\n", - " final_bits='',\n", - " final_wcs=True,\n", - " final_scale=scale)" + "Here we compare the single drizzled frames. Note that WF4 and WF2 overlap in the footprint on the sky." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "drz1_single = fits.getdata('u5jx010hr_single_sci.fits')\n", + "drz2_single = fits.getdata('u5jx0108r_single_sci.fits')\n", + "\n", + "fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(20, 10))\n", + "im1 = ax1.imshow(drz1_single, cmap='gray', vmin=0, vmax=4, origin='lower')\n", + "ax1.set_title('Target on WF4')\n", + "im2 = ax2.imshow(drz2_single, cmap='gray', vmin=0, vmax=4, origin='lower')\n", + "ax2.set_title('Target on WF2')\n", + "\n", + "x1 = ax1.get_position().get_points().flatten()[0]\n", + "x2 = ax2.get_position().get_points().flatten()[2] - x1\n", + "ax_cbar = fig.add_axes([x1, 0, x2, 0.03])\n", + "plt.colorbar(im1, cax=ax_cbar, orientation='horizontal')\n", + "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Display The Results of the Simulation Side-by-Side " + "### Inspect the single drizzled output frames derived from the dithered `c0m` inputs\n" ] }, { @@ -472,38 +817,179 @@ "metadata": {}, "outputs": [], "source": [ - "%matplotlib inline\n", - "drz_noneq = fits.getdata('nonequalized_sci.fits')\n", - "drz_eq = fits.getdata('equalized_sci.fits')\n", + "drz1_c0m_single = fits.getdata('u5jx010hr_c0m_single_sci.fits')\n", + "drz2_c0m_single = fits.getdata('u5jx0108r_c0m_single_sci.fits')\n", "\n", - "fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(24, 10))\n", - "im1 = ax1.imshow(drz_noneq, cmap='gray', vmin=0.88, vmax=1.08, origin='lower')\n", - "ax1.set_title('Original Simulation (Not equalized)')\n", - "im2 = ax2.imshow(drz_eq, cmap='gray', vmin=0.88, vmax=1.08, origin='lower')\n", - "ax2.set_title('Equalized Simulated Image')\n", + "fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(20, 10))\n", + "im1 = ax1.imshow(drz1_single, cmap='gray', vmin=0, vmax=4, origin='lower')\n", + "ax1.set_title('Target on WF4')\n", + "im2 = ax2.imshow(drz2_single, cmap='gray', vmin=0, vmax=4, origin='lower')\n", + "ax2.set_title('Target on WF2')\n", "\n", "x1 = ax1.get_position().get_points().flatten()[0]\n", "x2 = ax2.get_position().get_points().flatten()[2] - x1\n", "ax_cbar = fig.add_axes([x1, 0, x2, 0.03])\n", - "plt.colorbar(im1, cax=ax_cbar, orientation='horizontal')" + "plt.colorbar(im1, cax=ax_cbar, orientation='horizontal')\n", + "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The equalized image now have a ~constant background noise between all of the detectors with different sensitivities compared to the orignial simulation that has not been equalized." + "### Inspect the ratio of the two drizzled images to see changes due to photometric normalization" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ratio_flt_c0m = drz1_single/drz1_c0m_single\n", + "fig = plt.figure(figsize=(7, 7))\n", + "img_ratio = plt.imshow(ratio_flt_c0m, vmin=0.97, vmax=1.03, cmap='Greys_r', origin='lower')\n", + "plt.colorbar(img_ratio, orientation='vertical')\n", + "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ + "### Compare the combined drizzled image and the weight image. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "flt_drw_sci = fits.getdata('wfpc2_flt_drw_sci.fits')\n", + "flt_drw_wht = fits.getdata('wfpc2_flt_drw_wht.fits')\n", + "fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(20, 10))\n", + "im1 = ax1.imshow(flt_drw_sci, cmap='gray', vmin=0, vmax=0.5, origin='lower')\n", + "ax1.set_title('SCI')\n", + "im2 = ax2.imshow(flt_drw_wht, cmap='gray', vmin=0, vmax=50, origin='lower')\n", + "ax2.set_title('WHT')\n", + "x1 = ax1.get_position().get_points().flatten()[0]\n", + "x2 = ax2.get_position().get_points().flatten()[2] - x1\n", + "ax_cbar = fig.add_axes([x1, 0, x2, 0.03])\n", + "plt.colorbar(im1, cax=ax_cbar, orientation='horizontal')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note the Planetary Chip (PC) has ~4x WHT because the native pixel scale is ~2x smaller." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scale_ratio = 0.0996/0.0455\n", + "area_ratio = scale_ratio**2 \n", + "print('')\n", + "print('PIXEL SCALE RATIO: WF/PC = (0.0996/0.0455) = ', scale_ratio)\n", + "print('PIXEL AREA RATIO: WF/PC = (0.0996/0.0455)**2 = ', area_ratio)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we compare the two drizzled images, zooming in to highlight where the two chips overlap. Note that the one on the right has a slightly more uniform background across the field of view. Finally, we then identify three isolated stars for photometry. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "c0m_drw_sci = fits.getdata('wfpc2_c0m_drw_sci.fits')\n", + "flt_drw_sci = fits.getdata('wfpc2_flt_drw_sci.fits')\n", + "\n", + "im1 = c0m_drw_sci[700:1600, 700:1600]\n", + "im2 = flt_drw_sci[700:1600, 700:1600]\n", + "\n", + "fig, ax = plt.subplots(1, 2, figsize=(15, 10))\n", + "norm1 = ImageNormalize(im1, vmin=0.1, vmax=8, stretch=LogStretch())\n", + "ax[0].imshow(im1, norm=norm1, cmap='gray', origin='lower')\n", + "ax[1].imshow(im2, norm=norm1, cmap='gray', origin='lower')\n", + "ax[0].set_title('C0M drz_sci (zoom)', fontsize=20)\n", + "ax[1].set_title('FLT drz_sci (zoom)', fontsize=20)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot the ratio of the combined drizzled images from FLT and C0M" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ratio_drw = flt_drw_sci/c0m_drw_sci\n", + "fig = plt.figure(figsize=(7, 7))\n", + "img_ratio = plt.imshow(ratio_drw, vmin=0.97, vmax=1.03, cmap='Greys_r', origin='lower')\n", + "plt.colorbar(img_ratio, orientation='vertical')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Perform aperture photometry on three stars in the WF2, WF4 overlap region\n", + "\n", + "Here we perform aperture photometry in a radius of 5 pixels to show the change in count rates in each drizzle-combined image." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "positions = [(1180, 1381), (1388, 1190), (1232, 1210)] # Coords are (Y+1, X+1)\n", + "aperture = CircularAperture(positions, r=5)\n", + "annulus_aperture = CircularAnnulus(positions, r_in=10, r_out=15)\n", + "\n", + "print('Aperture radius:', aperture)\n", + "print('Annulus:', annulus_aperture)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "data = fits.getdata('wfpc2_flt_drw_sci.fits')\n", + "aperstats = ApertureStats(data, annulus_aperture)\n", + "bkg_mean = aperstats.mean\n", + "\n", + "phot_table = aperture_photometry(data, aperture)\n", "\n", - " Author: M. Cara, STScI Data Analysis Tools Branch\n", - " K. Huynh, WFC3 Instrument Branch\n", - " Created: December 14, 2018\n", - " Updated: November 16, 2023" + "for col in phot_table.colnames:\n", + " phot_table[col].info.format = '%.8g' # for consistent table output\n", + "\n", + "print('drz_flt_combined: \"wfpc2_flt_drw_sci.fits\"')\n", + "print(' SKY values :', bkg_mean)\n", + "print('---------------------------------')\n", + "print(phot_table)" ] }, { @@ -511,7 +997,83 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "data = fits.getdata('wfpc2_c0m_drw_sci.fits')\n", + "aperstats = ApertureStats(data, annulus_aperture)\n", + "bkg_mean = aperstats.mean\n", + "\n", + "phot_table = aperture_photometry(data, aperture)\n", + "\n", + "for col in phot_table.colnames:\n", + " phot_table[col].info.format = '%.8g' # for consistent table output\n", + "\n", + "print('drz_c0m_combined: \"wfpc2_c0m_drw_sci.fits\"')\n", + "print(' SKY values :', bkg_mean)\n", + "print('---------------------------------')\n", + "print(phot_table)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "While the total change in photometry is only a few DN/sec, properly accounting for the relative gain across chips is necessary to achieve the best possible photometric precision. As we've demonstrated in this notebook, this is accomplished by using the new `flt` style files, or running photeq on the older `c0m` style files." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Conclusions \n", + "[Table of Contents](#toc)\n", + "\n", + "MAST now include new WFPC2 `flt.fits` data products to replace the old `c0m.fits`, `c1m.fits` files. The new products include a new flux normalization between chips so that only a single zeropoint (inverse sensitivity, PHOTFLAM value) is needed for all chips. This allows for different chips to be drizzled together, for example, when combining images at multiple orientations. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " \n", + "## About this Notebook\n", + "[Table of Contents](#toc)\n", + "\n", + "***\n", + "\n", + " Created: 14 Dec 2018; M. Cara\n", + " Updated: 16 Nov 2023; K. Huynh & J. Mack\n", + " Updated: 10 Dec 2024; J. Mack & M. Revalski\n", + " \n", + " **Source:** [https://github.com/spacetelescope/hst_notebooks](https://github.com/spacetelescope/hst_notebooks)\n", + "\n", + "\n", + "### Additional Resources\n", + "\n", + "Below are some additional resources that may be helpful. Please send any questions to the [HST Help Desk](https://stsci.service-now.com/hst).\n", + "\n", + "- [WFPC2 Website](https://www.stsci.edu/hst/instrumentation/legacy/wfpc2)\n", + "- [WFPC2 2002 Data Handbook](https://www.stsci.edu/instruments/wfpc2/Wfpc2_dhb/WFPC2_longdhbcover.html) \n", + "- [WFPC2 2008 Instrument Handbook](https://www.stsci.edu/files/live/sites/www/files/home/hst/instrumentation/legacy/wfpc2/_documents/wfpc2_ihb.pdf)\n", + "- [Drizzlepac Handbook](https://hst-docs.stsci.edu/drizzpac)\n", + "- [WFPC2 Prior Drizzling Example](https://www.stsci.edu/files/live/sites/www/files/home/scientific-community/software/drizzlepac/examples/_documents/DrizzlePac_EX7.pdf)\n", + "\n", + "\n", + "\n", + "### Software Citations\n", + "If you use Python packages such as `astropy`, `astroquery`, `drizzlepac`, `matplotlib`, or `numpy` for published research, please cite the authors. Follow the links below for more information about citing various packages:\n", + "\n", + "* [Citing `astropy`](https://www.astropy.org/acknowledging.html)\n", + "* [Citing `astroquery`](https://github.com/astropy/astroquery/blob/main/astroquery/CITATION)\n", + "* [Citing `drizzlepac`](https://drizzlepac.readthedocs.io/en/latest/getting_started/citing_drizzlepac.html)\n", + "* [Citing `matplotlib`](https://matplotlib.org/stable/users/project/citing.html)\n", + "* [Citing `numpy`](https://numpy.org/citing-numpy/)\n", + "* [Citing `photutils`](https://photutils.readthedocs.io/en/stable/getting_started/citation.html)\n", + "***\n", + "\n", + "[Top of Page](#top)\n", + "\"Space" + ] } ], "metadata": { @@ -530,7 +1092,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.5" + "version": "3.11.7" }, "varInspector": { "cols": { diff --git a/notebooks/DrizzlePac/drizzle_wfpc2/requirements.txt b/notebooks/DrizzlePac/drizzle_wfpc2/requirements.txt index 64b6757f7..ea3c8adcb 100644 --- a/notebooks/DrizzlePac/drizzle_wfpc2/requirements.txt +++ b/notebooks/DrizzlePac/drizzle_wfpc2/requirements.txt @@ -1,12 +1,10 @@ -astropy>=5.3.3 -astroquery>=0.4.6 -drizzlepac>=3.5.1 -matplotlib>=3.7.0 -numpy>=1.23.4 -stsci.image>=2.3.5 -stsci.imagestats>=1.6.3 -stsci.skypac>=1.0.9 -stsci.stimage>=0.2.6 -stsci.tools>=4.0.1 -stwcs>=1.7.2 -crds +astropy>=6.0.1 +astroquery>=0.4.7 +crds>=11.17.20 +drizzlepac>=3.6.2 +ipython>=8.22.2 +jupyter>=1.0.0 +matplotlib>=3.8.4 +numpy>=1.26.4 +photutils>=1.12.0 +stwcs>=1.7.2 \ No newline at end of file