Skip to content

Commit

Permalink
Small edits so NB can be published Jan 2025.
Browse files Browse the repository at this point in the history
  • Loading branch information
@rjavila @gibsongreen
rjavila authored and gibsongreen committed Jan 8, 2025
1 parent 3ce6fe7 commit b07dc15
Showing 1 changed file with 91 additions and 37 deletions.
128 changes: 91 additions & 37 deletions notebooks/DrizzlePac/optimize_image_sampling/optimize_image_sampling.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,18 +21,14 @@
"[Import Packages](#imports)<br>\n",
"\n",
"[1. Retrieve Observations](#retrieve) <br>\n",
"\n",
"[2. Check image header data](#inspection) <br>\n",
"&nbsp;&nbsp;&nbsp;&nbsp;[2.1 Inspect useful keywords](#inspect_keywords)<br>\n",
"&nbsp;&nbsp;&nbsp;&nbsp;[2.2 Inspect the dither pattern](#inspect_dither)<br>\n",
"&nbsp;&nbsp;&nbsp;&nbsp;[2.3 Inspect the alignment](#inspect_wcs)<br>\n",
"\n",
"&nbsp;&nbsp;&nbsp;&nbsp;[2.3 Inspect the MAST alignment](#inspect_wcs)<br>\n",
"&nbsp;&nbsp;&nbsp;&nbsp;[2.4 Reset the WCS solution](#reset_wcs)<br>\n",
"[3. Drizzling](#drizzle) <br>\n",
"\n",
"[4. Results](#results) <br>\n",
"\n",
"[4. Inspecting the Results](#results) <br>\n",
"[5. Optimizing the final_pixfrac parameter](#pixfrac) <br>\n",
"\n",
"[6. Final Thoughts](#final) <br>\n",
"\n",
"[About this notebook](#about)<br>\n",
Expand Down Expand Up @@ -73,7 +69,8 @@
" - [**matplotlib**](https://matplotlib.org/stable/tutorials/pyplot.html) - make figures and graphics\n",
" - [**astropy**](https://www.astropy.org) - file handling, tables, units, WCS, statistics\n",
" - [**astroquery**](https://astroquery.readthedocs.io/en/latest/) - download data and query databases\n",
" - [**drizzlepac**](https://www.stsci.edu/scientific-community/software/drizzlepac) - align and combine HST images"
" - [**drizzlepac**](https://www.stsci.edu/scientific-community/software/drizzlepac) - align and combine HST images\n",
" - [**stwcs**] - tools for manipulating headerlets"
]
},
{
Expand All @@ -97,6 +94,7 @@
"from astropy import wcs\n",
"from astropy.table import Table\n",
"\n",
"from stwcs.wcsutil import headerlet\n",
"import drizzlepac\n",
"\n",
"%matplotlib inline\n",
Expand Down Expand Up @@ -200,6 +198,7 @@
"outputs": [],
"source": [
"flt_files = sorted(glob.glob('*fl?.fits'))\n",
"\n",
"data = []\n",
"keywords_ext0 = [\"ROOTNAME\", \"ASN_ID\", \"TARGNAME\", \"DETECTOR\", \"FILTER\", \"EXPTIME\", \n",
" \"RA_TARG\", \"DEC_TARG\", \"POSTARG1\", \"POSTARG2\", \"DATE-OBS\"]\n",
Expand Down Expand Up @@ -270,10 +269,9 @@
"metadata": {},
"source": [
"<a id=\"inspect_wcs\"></a>\n",
"### 2.3 Inspect the alignment \n",
"### 2.3 Inspect the MAST alignment \n",
"\n",
"Check the active WCS solution in the image header. If the image is aligned to a catalog, list the number of matches and the fit RMS in mas. <br> \n",
"Convert the fit RMS values to pixels for comparison with the alignment results performed later in this notebook."
"Check the active WCS solution in the FITS image header. If the image is aligned to a catalog, list the number of matches and the fit RMS in milli-arcseconds and in pixels."
]
},
{
Expand Down Expand Up @@ -318,9 +316,62 @@
"source": [
"Here we see that the four FLC images have a \"FIT_REL_GSC242\" WCS solution, which means they were aligned as a set to the reference catalog 'GSC v2.4.2' using the drizzled-combined (DRC) image, and that WCS was propogated back to the original FLCs. This is reflected in the fact that the number of matches and fit rms values are the same for each FLC frame. \n",
"\n",
"For more details on absolute astrometry in HST images, see [Section 4.5 in the DrizzlePac Handbook](https://hst-docs.stsci.edu/drizzpac/chapter-4-astrometric-information-in-the-header/4-5-absolute-astrometry)\n",
"Note that the alignment was based on 26 matches to the GSC v2.4.2 catalog. This means that the target did not successfully align to the Gaia eDR3 solution since there are very few stars in the frame. \n",
"\n",
"We can various functions within the `headerlet` of the `stwcs` package to find, list, and ultimately update the WCS solutions in an image. Let's make a list of the images and then print out all the WCS solutions in the first image in the list."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"wcsnames = headerlet.get_headerlet_kw_names(flt_files[0], kw='WCSNAME')\n",
"wcsnames"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For more details on the WCS naming convention, see [DrizzlePac readthedocs](https://drizzlepac.readthedocs.io/en/latest/mast_data_products/astrometry.html#interpreting-wcs-names)\n",
"\n",
"For more discussion on absolute astrometry, see [Section 4.5 in the DrizzlePac Handbook](https://hst-docs.stsci.edu/drizzpac/chapter-4-astrometric-information-in-the-header/4-5-absolute-astrometry)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id=\"reset_wcs\"></a>\n",
"### 2.4 Reset the WCS solution\n",
"\n",
"Next, we find the index of the WCS solution we want to use, and then apply that solution to each of the science images. Because the program was designed to have precise dithers, we reset the WCS to the “a priori” (distortion-only) solution to preserve relative alignment and get the sharpest PSFs when drizzling the images together.\n",
"\n",
"In this example, the \"a priori\" solution is named **IDC_w3m18525i-GSC240**, in which the coordinates of the guide stars were corrected to the coordinates reported in GAIA DR1, with no additional catalog fitting.\n",
"\n",
"Since the relative alignment of the individual FLC frames is preserved, we do not need run `TweakReg` before drizzling the images."
"Here we first grab the index of list element with value 'IDC_w3m18525i-GSC240'. The index in this list + 1 is the same as the EXTVER of the corresponding HDRLET. We need to add 1 because lists are 0-indexed, while EXTVER's are 1 indexed."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for flt_file in flt_files:\n",
" \n",
" chosen_ext = wcsnames.index('IDC_w3m18525i-GSC240')+1\n",
" headerlet.restore_from_headerlet(flt_file, hdrext=('HDRLET', chosen_ext), archive=False, force=False)\n",
" print(f'Image: {flt_file} WCSNAME: {fits.getval(flt_file,\"wcsname\",ext=1)}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since the relative alignment of the individual FLC frames is preserved with the 'a priori' solution, we do not need run `TweakReg` before drizzling the images. "
]
},
{
Expand All @@ -331,10 +382,12 @@
"## 3. Drizzling \n",
"[Table of Contents](#toc)\n",
"\n",
"Before drizzling can be performed, a plate scale for the output image should be chosen. This will usually be dictated by the science needs. In theory, critical sampling of a PSF occurs with 2.355 pixels. The FWHM of the WFC3/IR detector at 1600 nm is ~0.151\" (Table 7.5 of [WFC Instrument Handbook](https://hst-docs.stsci.edu/wfc3ihb/chapter-7-ir-imaging-with-wfc3/7-6-ir-optical-performance)). Using these values, a plate scale of 0.065\"/pix will be selected for the drizzled output frame, approximately half of the native plate scale (0.1283\"/pix). \n",
"Before drizzling can be performed, a plate scale for the output image should be chosen. This will usually be dictated by the science needs. In theory, critical sampling of a PSF occurs with 2.355 pixels. The FWHM of the WFC3/IR detector at 1600 nm is ~0.151\" (Table 7.5 of [WFC Instrument Handbook](https://hst-docs.stsci.edu/wfc3ihb/chapter-7-ir-imaging-with-wfc3/7-6-ir-optical-performance)). \n",
"\n",
"<b>Using these values, we select a plate scale of 0.065\"/pix for the drizzled output frame, approximately half of the native plate scale (0.1283\"/pix). \n",
"\n",
"The mechanism by which the *drizzle* algorithm improves sampling in the output image is by shrinking the input pixels before adding them to the output pixel grid (see Figure 2 in [Fruchter and Hook (2002)](http://iopscience.iop.org/article/10.1086/338393/pdf)). In practice this process is controlled in the `AstroDrizzle` package by adjusting the `final_pixfrac` parameter, which is the fractional size of the pixel to be used. \n",
"Below `AstroDrizzle` is called twice. The first time, no changes are made to the plate scale or the size of the drop. In the second call, the plate scale and pixfrac are reduced. For a four point dither stack, a `final_pixfrac` value of 0.8 is usually optimal.\n",
"Below `AstroDrizzle` is called twice. The first time, no changes are made to the plate scale or the size of the drop. In the second call, the plate scale and pixfrac are reduced. \n",
"\n",
"Next, we get some recommended values for drizzling from the MDRIZTAB reference file. The parameters in this file are different for each detector and are based on the number of input frames. These are a good starting point for drizzling and may be adjusted accordingly.\n",
"\n",
Expand Down Expand Up @@ -430,7 +483,7 @@
"metadata": {},
"source": [
"<a id=\"results\"></a>\n",
"## 4. Results \n",
"## 4. Inspecting the Results \n",
"[Table of Contents](#toc)\n",
"\n",
"The drizzled science and weight images produced from the first call to AstroDrizzle with no optimization of the plate scale and pixfrac are plotted below. "
Expand Down Expand Up @@ -515,7 +568,7 @@
"\n",
"Below, a series of experiments will be run to determine the best `final_pixfrac` value for the selected output plate scale, by varying pixfrac in steps of 0.1 over a range of values from 0.1 to 1.0.\n",
"\n",
"<div class=\"alert alert-block alert-warning\" style=\"color:black\" > This cell may take a few minutes to execute. </div>"
"<div class=\"alert alert-block alert-warning\" style=\"color:black\" > This cell may take a several minutes to execute. </div>"
]
},
{
Expand All @@ -533,9 +586,11 @@
" drizzlepac.astrodrizzle.AstroDrizzle(flt_files,\n",
" output=outname,\n",
" **selected_params,\n",
" context=False, build=True, preserve=False, \n",
" context=False, \n",
" build=True, \n",
" preserve=False, \n",
" skymethod='match',\n",
" driz_separate=False, median=False, blot=False, driz_cr=False,\n",
" driz_separate=False, median=False, blot=False, driz_cr=False, # CR-rej off\n",
" final_pixfrac=pixfrac,\n",
" final_wcs=True,\n",
" final_rot=0.,\n",
Expand All @@ -562,7 +617,7 @@
" sci1 = hdu1[1].data\n",
" wht1 = hdu1[2].data\n",
" \n",
"with fits.open('f160w_0.8_drz.fits') as hdu2:\n",
"with fits.open('f160w_0.6_drz.fits') as hdu2:\n",
" sci2 = hdu2[1].data\n",
" wht2 = hdu2[2].data\n",
" \n",
Expand All @@ -577,7 +632,7 @@
"ax[0, 0].imshow(sci1, norm=norm2, cmap='gray')\n",
"ax[0, 0].text(1640, 1920, 'final_pixfrac=0.1', fontsize='20', color='maroon')\n",
"ax[0, 1].imshow(sci2, norm=norm2, cmap='gray')\n",
"ax[0, 1].text(1640, 1920, 'final_pixfrac=0.8', fontsize='20', color='maroon')\n",
"ax[0, 1].text(1640, 1920, 'final_pixfrac=0.6', fontsize='20', color='maroon')\n",
"ax[0, 2].imshow(sci3, norm=norm2, cmap='gray')\n",
"ax[0, 2].text(1640, 1920, 'final_pixfrac=1.0', fontsize='20', color='maroon')\n",
"ax[1, 0].imshow(wht1, norm=norm3, cmap='gray')\n",
Expand All @@ -591,7 +646,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"This figure above shows the central region of the science and weight images produced by `final_pixfrac` values (left to right) of 0.1, 0.8 and 1.0. The top row is the science frame, the bottom is the weight image. The science image with the smallest `final_pixfrac` value of 0.1 shows a noisy background and holes in the image where no input pixels fall into the output grid because pixfrac is too small. Inspection of the weight map corresponding to the smallest `final_pixfrac` value shows many places with weights of zero, indicating that a `final_pixfrac` value of 0.1 is clearly too small and was only included in this example for illustrative purposes."
"This figure above shows the central region of the science and weight images produced by `final_pixfrac` values (left to right) of 0.1, 0.8 and 1.0. The top row is the science frame, the bottom is the weight image. The science image with the smallest `final_pixfrac` value of 0.1 shows a noisy background and holes in the image where no input pixels fall into the output grid because pixfrac is too small. \n",
"\n",
"Inspection of the weight map corresponding to the smallest `final_pixfrac` value shows many places with weights of zero, indicating that a `final_pixfrac` value of 0.1 is clearly too small and was only included in this example for illustrative purposes."
]
},
{
Expand All @@ -600,6 +657,8 @@
"source": [
"Another piece of information that can be useful is the amount of noise in the weight image. As suggested in the [DrizzlePac Handbook](https://hst-docs.stsci.edu/drizzpac/chapter-7-data-quality-checks-and-trouble-shooting-problems/7-3-inspecting-drizzled-products-after-user-reprocessing) Section 7.3, statistics performed on the drizzled weight image should yield a RMS/median value less than ~0.2. This threshold controls the trade-off between improving image resolution versus increasing background noise due to pixel resampling.\n",
" \n",
"The final_pixfrac value should be small enough to avoid degrading the final drizzle-combined image, but large enough that when all images are \"dropped\" onto the final frame, coverage of the output frame is fairly uniform. In general, final_pixfrac should be slightly larger than the final output scale to allow some \"spillover\" to adjacent pixels. This will help avoid \"holes\" in the final product when a given pixel has been flagged as bad in several frames.\n",
" \n",
"The figure below shows the RMS/median as a function of `final_pixfrac`. One should take care to use the same region in the weight image as the region where the object of interest is located in the science image. If one is using the entire image for scientific analysis, then one should measure the statistics of the weight image where there is more variance. "
]
},
Expand Down Expand Up @@ -633,21 +692,23 @@
"# Plotting commands \n",
"plt.clf()\n",
"plt.xlim(1.025, -0.025)\n",
"plt.ylim(0., 0.9)\n",
"plt.ylim(0., 0.3)\n",
"plt.scatter(fraclist, std_med, s=50)\n",
"plt.axhline(0.2, ls='--', lw=3, c='r')\n",
"plt.axhline(0.2, ls='--', lw=3, c='g')\n",
"plt.xlabel('pixfrac', fontsize=18)\n",
"plt.ylabel('rms/median', fontsize=18)\n",
"plt.text(1.0, 0.85, target, fontsize=16, horizontalalignment='left')\n",
"plt.text(1.0, 0.80, 'Scale=' + str(scale) + '\"', fontsize=16, horizontalalignment='left')\n",
"plt.text(1.0, 0.75, str(nimg) + ' images', fontsize=16, horizontalalignment='left')"
"plt.text(0.95, 0.28, target, fontsize=14, horizontalalignment='left')\n",
"plt.text(0.95, 0.26, 'Scale=' + str(scale) + '\"', fontsize=14, horizontalalignment='left')\n",
"plt.text(0.95, 0.24, str(nimg) + ' images', fontsize=14, horizontalalignment='left')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The RMS/median increases steadily down to `final_pixfrac` value of 0.3, where there is a change in the slope of the function, indicating that something odd might be happening beyond that point. The RMS remains below 20% of the median all the way down to `final_pixfrac` of 0.8, but the science and weight images should be inspected at all values to ensure no quality degradation has occurred. "
"Note that none of the test images have RMS/median greater than 0.2, so this does not seem to be a good metric for WFC3/IR in which the PSF is undersampled. (We suspect that this recommendation from the DrizzlePac Handbook was based on WFC3/UVIS or ACS/WFC dithered frames.)\n",
"\n",
"The RMS/median increases steadily to a `final_pixfrac` value of 0.5, where there is a change in the slope of the function. For an output scale of 0.065\"/pix (0.5 times the native plate scale) thus select the pixfrac 0.6 for the final drizzled image. "
]
},
{
Expand All @@ -663,7 +724,7 @@
"## 6. Final thoughts \n",
"[Table of Contents](#toc)\n",
"\n",
"The `final_pixfrac` value has to be small enough to avoid degrading the combined output image, but large enough that when all images are \"dropped\" onto the final frame, coverage of the output frame is fairly uniform. In general, `final_pixfrac` should be slightly larger than the final output scale to allow some \"spillover\" to adjacent pixels. This will help avoid \"holes\" in the final product when a given pixel has been flagged as \"bad\" in several frames. As a rule of thumb, statistics performed on the drizzled weight image in the region of interest should yield an RMS value (standard deviation) that is less than 20% of the median value. This threshold is a balance between the benefits of improving the image resolution at the expense of increasing noise in the background."
"The `final_pixfrac` value has to be small enough to avoid degrading the combined output image, but large enough that when all images are \"dropped\" onto the final frame, coverage of the output frame is fairly uniform. In general, `final_pixfrac` should be slightly larger than the final output scale to allow some \"spillover\" to adjacent pixels. This will help avoid \"holes\" in the final product when a given pixel has been flagged as \"bad\" in several frames. Optimizing the final drizzle parameters is a balance between the benefits of improving the image resolution and PSFs at the expense of increasing noise in the background."
]
},
{
Expand All @@ -674,7 +735,7 @@
"## About this Notebook\n",
"\n",
" Created: 12 Dec 2018; R. Avila\n",
" Updated: 29 Apr 2024; R. Avila & J. Mack\n",
" Updated: 30 Dec 2024; R. Avila & J. Mack\n",
"\n",
"**Source:** GitHub [spacetelescope/hst_notebooks](https://github.com/spacetelescope/hst_notebooks)\n",
"\n",
Expand Down Expand Up @@ -710,13 +771,6 @@
"source": [
"[Top of Page](#top)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand Down

0 comments on commit b07dc15

Please sign in to comment.