<jupyter_start><jupyter_text>Aligning HST images to an absolute reference catalog*** Note: The notebook in this repository 'Initializtion.ipynb' goes over many of the basic concepts such as the setup of the environment/package installation and should be read first if you are new to HST images, 'DrizzlePac' or 'Astroquery'. Note: This notebook is based on WFC3 ISR 2017-19: Aligning HST Images to Gaia: a Faster Mosaicking Workflow and contains a subset of the information/code found in the repository here.  For more information, see the notebook in that repository titled 'Gaia_alignment.ipynb'.  Introduction The alignment of HST exposures is a critical step in image stacking/combination performed by software such as `AstroDrizzle`.  Generally, a relative alignment is performed that aligns one image (or multiple images) to another image which is designated as the reference image.  This makes it so the images are aligned to each other, but the pointing error of the observatory can still cause the images to have incorrect absolute astrometry.When absolute astrometry is desired, the images can be aligned to an external catalog that is known to be on an absolute frame.  In this example, we will provide a workflow to query catalogs such as SDSS and Gaia via the astroquery package, and then align the images to that catalog via TweakReg.For more information about TweakReg, see the other notebooks in this repository or the __[TweakReg Documentation](https://drizzlepac.readthedocs.io/en/deployment/tweakreg.html)__.For more information on Astroquery, see the other notebooks in this repository or the __[Astroquery Documentation](https://astroquery.readthedocs.io/en/latest/)__.<jupyter_code>import astropy.units as u
import glob
import numpy as np
import matplotlib.pyplot as plt
import os

from astropy.io import fits
from astropy.table import Table
from astropy.units import Quantity
from astropy.coordinates import SkyCoord
from astroquery.gaia import Gaia
from astroquery.mast import Observations
from astroquery.sdss import SDSS

from ccdproc import ImageFileCollection
from IPython.display import Image

from drizzlepac import tweakreg
from drizzlepac import astrodrizzle<jupyter_output><empty_output><jupyter_text>1. Download the DataFor this example, we will use WFC3/UVIS images of NGC 6791 from Visit 01 of proposal 12692.<jupyter_code># Get the observation records
obsTable = Observations.query_criteria(obs_id='ibwb01*', proposal_id=12692, obstype='all', filters='F606W')

# Get the listing of data products
products = Observations.get_product_list(obsTable)

# Filter the products for exposures
filtered_products = Observations.filter_products(products, productSubGroupDescription='FLC')

# Show the table
filtered_products
# Download all the images above
Observations.download_products(filtered_products, mrp_only=False)
# For convenience, move the products into the current directory.
for flc in glob.glob('./mastDownload/HST/*/*flc.fits'):
    flc_name = os.path.split(flc)[-1]
    os.rename(flc, flc_name)<jupyter_output><empty_output><jupyter_text>Inspect the image headerThe cell below shows how to query information from the image header using `ImageFileCollection` in `ccdproc`. We see that the 1st exposure is 30 seconds and the 2nd and 3rd exposures are 360 seconds. The 3rd exposure is dithered by ~82" in the Y-direction which is approximately the width of one UVIS chip.<jupyter_code>collec = ImageFileCollection('./', glob_include="*flc.fits", ext=0,
                             keywords=["targname", "ra_targ", "dec_targ", "filter", "exptime", "postarg1", "postarg2"])

table = collec.summary
table['exptime'].format = '7.1f'
table['ra_targ'].format = '7.7f'
table['dec_targ'].format = '7.7f'
table['postarg1'].format = '7.2f'
table['postarg2'].format = '7.2f'
table<jupyter_output><empty_output><jupyter_text>2. Querying catalogsNow that we have the images, we will download the reference catalogs from both SDSS and Gaia using `astroquery`.  2a. Identify CoordinatesWe will first create a SkyCoord Object to point astroquery to where we are looking on the sky.  Since our example uses data from NGC 6791, we will use the `ra_targ` and `dec_targ` keywords from the first image to get the coordinates of the object.<jupyter_code>RA = table['ra_targ'][0]
Dec = table['dec_targ'][0]<jupyter_output><empty_output><jupyter_text>2b. SDSS QueryWe now give those values to an astropy `SkyCoord` object, which we will pass to the SDSS.  Additionally, we use an astropy `Quantity` object to create a radius for the SDSS query.  We set the radius to 6 arcminutes to comfortably cover the area of our images. For reference UVIS detector field of view is ~2.7'x2.7' and a y-dither of 82" covers a total area on the sky of ~2.7'x4.1'.<jupyter_code>coord = SkyCoord(ra=RA, dec=Dec, unit=(u.deg, u.deg))
radius = Quantity(6., u.arcmin)<jupyter_output><empty_output><jupyter_text>Then we only need to perform the query via the `SDSS.query_region` method of `astroquery.sdss`. The `spectro=False` keyword argument means we want to exclude spectroscopic objects, as we are looking for objects to match with an image.  In the fields parameter, we specify a list of fields we want returned by the query.  In this case we only need the position, and maybe a magnitude 'g' if we want to cut very dim and/or bright objects out of the catalog, as those are likely measured poorly. Details on selecting objects by magnitude may be found in the original ['Gaia_alignment' notebook](https://github.com/spacetelescope/gaia_alignment).  Many other fields are available in the SDSS query and are [documented here](http://cas.sdss.org/dr7/en/help/browser/description.asp?n=PhotoObj&t=V).<jupyter_code>sdss_query = SDSS.query_region(coord, radius=radius, spectro=False, fields=['ra', 'dec', 'g'])
sdss_query<jupyter_output><empty_output><jupyter_text>We then just need to save the catalog to use with TweakReg.  Since the returned value of the query is an `astropy.Table`, saving the file is very straightforward:<jupyter_code>sdss_query.write('sdss.cat', format='ascii.commented_header')<jupyter_output><empty_output><jupyter_text>2c. Gaia QuerySimilarly to SDSS, we can query Gaia catalogs for our target via `astroquery.gaia`.  We can use the same `coord` and  `radius` from the SDSS query.<jupyter_code>gaia_query = Gaia.query_object_async(coordinate=coord, radius=radius)
gaia_query<jupyter_output><empty_output><jupyter_text>This query has returned very large number of columns. We want to pare down the catalog to make it easier to use with `TweakReg`.  We can select only the useful columns via:<jupyter_code>reduced_query = gaia_query['ra', 'dec', 'phot_g_mean_mag']
reduced_query<jupyter_output><empty_output><jupyter_text>Then we write this catalog to an ascii file for use with `TweakReg`.<jupyter_code>reduced_query.write('gaia.cat', format='ascii.commented_header')<jupyter_output><empty_output><jupyter_text>3. Running TweakRegWith the catalogs downloaded and the headers populated, we simply need to run TweakReg with each catalog passed into the `refcat` parameter. The steps below compare the astrometric residuals obtained from aligning to each `refcat`. In each test case, we set `updatehdr` to False until we are satisfied with the alignment by inspecting both the shift file and the astrometric residual plots. 3a. SDSS Alignment<jupyter_code>refcat = 'sdss.cat'
cw = 3.5  # Set to two times the FWHM of the PSF of the UVIS detector
wcsname = 'SDSS'  # Specify the WCS name for this alignment

tweakreg.TweakReg('*flc.fits',  # Pass input images
                  updatehdr=False,  # update header with new WCS solution
                  imagefindcfg={'threshold': 500., 'conv_width': cw},  # Detection parameters, threshold varies for different data
                  refcat=refcat,  # Use user supplied catalog (Gaia)
                  interactive=False,
                  see2dplot=False,
                  shiftfile=True,  # Save out shift file (so we can look at shifts later)
                  outshifts='SDSS_shifts.txt',  # name of the shift file
                  wcsname=wcsname,  # Give our WCS a new name
                  reusename=True,
                  ylimit=1.5,
                  fitgeometry='general')  # Use the 6 parameter fit<jupyter_output><empty_output><jupyter_text>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:- Filename- X Shift [pixels]- Y Shift [pixels]- Rotation [degrees]- Scale- X RMS [pixels]- Y RMS [pixels]<jupyter_code>for line in open('SDSS_shifts.txt').readlines():
    print(line)
# Astrometric residual plots
Image(filename='residuals_ibwb01xqq_flc.png',width=500, height=300)
Image(filename='residuals_ibwb01xrq_flc.png',width=500, height=300)
Image(filename='residuals_ibwb01xxq_flc.png',width=500, height=300)<jupyter_output><empty_output><jupyter_text>As we can see, the RMS is fairly large at about 0.5 pixels, which is not a great fit. This is likely because the SDSS astrometric precision is not high enough to get good HST alignment. One approach would be to align the first image to SDSS and then align the remaining HST images to one another. This would improve both the absolute and relative alignment of the individual frames. 3b. Gaia Alignment<jupyter_code>refcat = 'gaia.cat'
cw = 3.5  # Set to two times the FWHM of the PSF.
wcsname = 'Gaia'  # Specify the WCS name for this alignment

tweakreg.TweakReg('*flc.fits',  # Pass input images
                  updatehdr=False,  # update header with new WCS solution
                  imagefindcfg={'threshold':500.,'conv_width':cw},  # Detection parameters, threshold varies for different data
                  refcat=refcat,  # Use user supplied catalog (Gaia)
                  interactive=False,
                  see2dplot=False,
                  shiftfile=True,  # Save out shift file (so we can look at shifts later)
                  outshifts='Gaia_shifts.txt',  # name of the shift file
                  wcsname=wcsname,  # Give our WCS a new name
                  reusename=True,
                  sigma=2.3,
                  ylimit=0.2,
                  fitgeometry='general')  # Use the 6 parameter fit<jupyter_output><empty_output><jupyter_text>We can similarly look at the shift file from alignment to the Gaia catalog:<jupyter_code>for line in open('Gaia_shifts.txt').readlines():
    print(line)
# Astrometric residual plots
Image(filename='residuals_ibwb01xqq_flc.png',width=500, height=300)
Image(filename='residuals_ibwb01xrq_flc.png',width=500, height=300)
Image(filename='residuals_ibwb01xxq_flc.png',width=500, height=300)<jupyter_output><empty_output><jupyter_text>As expected, the Gaia catalog does quite a bit better, with rms residuals less tha 0.05 pixels.  To apply these transformations to the image, we simply need to run TweakReg the same as before, but set the parameter `updatehdr` equal to `True`:<jupyter_code>refcat = 'gaia.cat'
cw = 3.5  # Set to two times the FWHM of the PSF.
wcsname = 'Gaia'  # Specify the WCS name for this alignment

tweakreg.TweakReg('*flc.fits',  # Pass input images
                  updatehdr=True,  # update header with new WCS solution
                  imagefindcfg={'threshold': 500., 'conv_width': cw},  # Detection parameters, threshold varies for different data
                  refcat=refcat,  # Use user supplied catalog (Gaia)
                  interactive=False,
                  see2dplot=False,
                  shiftfile=True,  # Save out shift file (so we can look at shifts later)
                  outshifts='Gaia_shifts.txt',  # name of the shift file
                  wcsname=wcsname,  # Give our WCS a new name
                  reusename=True,
                  sigma=2.3,
                  fitgeometry='general')  # Use the 6 parameter fit<jupyter_output><empty_output><jupyter_text>4. Drizzle the DataWhile the three sets of FLC files are now aligned, we drizzle together only the two long exposures. When exposures are very different lengths, drizzling them together doesn't work well when 'EXP' weighting is used. For objects that saturate in the long exposures, the problem occurs at the boundary where you transition from only short exposure to short plus long. Here the pixels getting power from long exposure pixels are only getting power from pixels whose centers are outside the ring, and thus they are weighted lower than they would be if they were getting values from both inside and outside the ring. The result is a discontinuity in the PSF radial profile and a resulting flux which is too low in those boundary pixels. For photometry of saturated objects, the short exposures should be drizzled separately from the long exposures.<jupyter_code>astrodrizzle.AstroDrizzle('ibwb01x[rx]q_flc.fits', 
                          output='f606w',
                          preserve=False,
                          clean=True, 
                          build=False,
                          context=False,
                          skymethod='match',
                          driz_sep_bits='64, 32',
                          combine_type='minmed',
                          final_bits='64, 32')
# Display the combined science and weight images 

sci = fits.getdata('f606w_drc_sci.fits')
wht = fits.getdata('f606w_drc_wht.fits')

fig = plt.figure(figsize=(20, 20))
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)

ax1.imshow(sci, vmin=-0.05, vmax=0.4, cmap='Greys_r', origin='lower')
ax2.imshow(wht, vmin=0, vmax=1000, cmap='Greys_r', origin='lower')<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup