.. only:: html
.. note::
:class: sphx-glr-download-link-note
Click :ref:`here <sphx_glr_download_examples_plot_c_miri.py>`
to download the full example code
.. rst-class:: sphx-glr-example-title
Aligning JWST/MIRI images with JHAT.
An example MIRI Dataset is downloaded, and then a series of alignment methods are used. For more information on the key parameters used for alignment see :ref:`params:Useful Parameters`.
import sys,os,glob
from astropy.io import fits
from astropy.table import Table
from astropy.nddata import extract_array
from astropy.coordinates import SkyCoord
from astropy import wcs
from astropy.wcs.utils import skycoord_to_pixel
from astropy import units as u
import numpy as np
import matplotlib.pyplot as plt
from astroquery.mast import Observations
from astropy.visualization import (simple_norm,LinearStretch)
import jhat
from jhat import jwst_photclass,st_wcs_align
Download some Data
For this example we download 2 MIRI cal images from MAST. They're the same field and different filters. Note that the code will also work for level 3 images.
obs_table1 = Observations.query_criteria(obs_id='jw02107-o038_t019_miri_f770w')
data_products_by_obs = Observations.get_product_list(obs_table1)
data_products_by_obs = data_products_by_obs[data_products_by_obs['calib_level']==2]
data_products_by_obs = data_products_by_obs[data_products_by_obs['productSubGroupDescription']=='CAL'][0]
Observations.download_products(data_products_by_obs,extension='fits')
obs_table2 = Observations.query_criteria(obs_id='jw02107-c1018_t019_miri_f1130w')
data_products_by_obs = Observations.get_product_list(obs_table2)
data_products_by_obs = data_products_by_obs[data_products_by_obs['calib_level']==2]
data_products_by_obs = data_products_by_obs[data_products_by_obs['productSubGroupDescription']=='CAL'][0]
Observations.download_products(data_products_by_obs,extension='fits')
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:JWST/product/jw02107038001_02101_00001_mirimage_cal.fits to ./mastDownload/JWST/jw02107038001_02101_00001_mirimage/jw02107038001_02101_00001_mirimage_cal.fits ... [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:JWST/product/jw02107038001_02105_00001_mirimage_cal.fits to ./mastDownload/JWST/jw02107038001_02105_00001_mirimage/jw02107038001_02105_00001_mirimage_cal.fits ... [Done]
Table length=1
| Local Path | Status | Message | URL |
|---|---|---|---|
| str98 | str8 | object | object |
| ./mastDownload/JWST/jw02107038001_02105_00001_mirimage/jw02107038001_02105_00001_mirimage_cal.fits | COMPLETE | None | None |
Examine the Reference Image
files = glob.glob('mastDownload/JWST/*miri*/*cal.fits')
ref_image = files[0]
print(ref_image)
ref_fits = fits.open(ref_image)
ref_data = fits.open(ref_image)['SCI',1].data
norm1 = simple_norm(ref_data,stretch='log',min_cut=5,max_cut=25)
plt.imshow(ref_data, origin='lower',
norm=norm1,cmap='gray')
plt.gca().tick_params(labelcolor='none',axis='both',color='none')
plt.show()
.. image-sg:: /examples/images/sphx_glr_plot_c_miri_001.png :alt: plot c miri :srcset: /examples/images/sphx_glr_plot_c_miri_001.png :class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
.. code-block:: none
mastDownload/JWST/jw02107038001_02101_00001_mirimage/jw02107038001_02101_00001_mirimage_cal.fits
Zoom in to see the offset
Here add an artificial offset to the wcs, and then we see the same star in both images at the same ra/dec location, demonstrating a large offset between the images.
star_location = SkyCoord('23:09:44.0809','-43:26:05.613',unit=(u.hourangle,u.deg))
align_image = files[1]
align_fits = fits.open(align_image)
align_fits['SCI',1].header['CRPIX1']+=2
align_fits['SCI',1].header['CRPIX2']+=2
align_fits.writeto(align_image,overwrite=True)
align_data = fits.open(align_image)['SCI',1].data
ref_y,ref_x = skycoord_to_pixel(star_location,wcs.WCS(ref_fits['SCI',1],ref_fits))
align_y,align_x = skycoord_to_pixel(star_location,wcs.WCS(align_fits['SCI',1],align_fits))
ref_cutout = extract_array(ref_data,(11,11),(ref_x,ref_y))
align_cutout = extract_array(align_data,(11,11),(align_x,align_y))
norm1 = simple_norm(ref_cutout,stretch='log',min_cut=-1,max_cut=200)
norm2 = simple_norm(align_cutout,stretch='log',min_cut=-1,max_cut=200)
fig,axes = plt.subplots(1,2)
axes[0].imshow(ref_cutout, origin='lower',
norm=norm1,cmap='gray')
axes[1].imshow(align_cutout, origin='lower',
norm=norm2,cmap='gray')
axes[0].set_title('Reference')
axes[1].set_title('To Align')
axes[0].tick_params(labelcolor='none',axis='both',color='none')
axes[1].tick_params(labelcolor='none',axis='both',color='none')
plt.show()
.. image-sg:: /examples/images/sphx_glr_plot_c_miri_002.png :alt: Reference, To Align :srcset: /examples/images/sphx_glr_plot_c_miri_002.png :class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
.. code-block:: none
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/astropy/wcs/wcs.py:725: FITSFixedWarning: 'datfix' made the change 'Set DATE-BEG to '2022-07-06T17:29:42.548' from MJD-BEG.
Set DATE-AVG to '2022-07-06T17:29:53.648' from MJD-AVG.
Set DATE-END to '2022-07-06T17:30:04.748' from MJD-END'.
warnings.warn(
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/astropy/wcs/wcs.py:725: FITSFixedWarning: 'obsfix' made the change 'Set OBSGEO-L to -72.176807 from OBSGEO-[XYZ].
Set OBSGEO-B to -38.353152 from OBSGEO-[XYZ].
Set OBSGEO-H to 1740801417.596 from OBSGEO-[XYZ]'.
warnings.warn(
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/astropy/wcs/wcs.py:725: FITSFixedWarning: 'datfix' made the change 'Set DATE-BEG to '2022-07-06T17:47:53.158' from MJD-BEG.
Set DATE-AVG to '2022-07-06T17:48:32.008' from MJD-AVG.
Set DATE-END to '2022-07-06T17:49:10.859' from MJD-END'.
warnings.warn(
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/astropy/wcs/wcs.py:725: FITSFixedWarning: 'obsfix' made the change 'Set OBSGEO-L to -72.174733 from OBSGEO-[XYZ].
Set OBSGEO-B to -38.353284 from OBSGEO-[XYZ].
Set OBSGEO-H to 1740817774.322 from OBSGEO-[XYZ]'.
warnings.warn(
Create a Photometric Catalog for Relative Alignment
We choose one of the images to be the reference image, and then create a catalog that we will use to align the other image.
jwst_phot = jwst_photclass()
jwst_phot.run_phot(imagename=ref_image,photfilename='auto',overwrite=True)
ref_catname = ref_image.replace('.fits','.phot.txt') # the default
refcat = Table.read(ref_catname,format='ascii')
print(refcat)
.. rst-class:: sphx-glr-script-out
.. code-block:: none
0 mastDownload/JWST/jw02107038001_02101_00001_mirimage/jw02107038001_02101_00001_mirimage_cal.phot.txt
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/astropy/wcs/wcs.py:725: FITSFixedWarning: 'datfix' made the change 'Set DATE-BEG to '2022-07-06T17:29:42.548' from MJD-BEG.
Set DATE-AVG to '2022-07-06T17:29:53.648' from MJD-AVG.
Set DATE-END to '2022-07-06T17:30:04.748' from MJD-END'.
warnings.warn(
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/astropy/wcs/wcs.py:725: FITSFixedWarning: 'obsfix' made the change 'Set OBSGEO-L to -72.176807 from OBSGEO-[XYZ].
Set OBSGEO-B to -38.353152 from OBSGEO-[XYZ].
Set OBSGEO-H to 1740801417.596 from OBSGEO-[XYZ]'.
warnings.warn(
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/astropy/stats/sigma_clipping.py:411: AstropyUserWarning: Input data contains invalid values (NaNs or infs), which were automatically clipped.
warnings.warn('Input data contains invalid values (NaNs or '
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/astropy/stats/sigma_clipping.py:411: AstropyUserWarning: Input data contains invalid values (NaNs or infs), which were automatically clipped.
warnings.warn('Input data contains invalid values (NaNs or '
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/astropy/units/function/logarithmic.py:47: RuntimeWarning: invalid value encountered in log10
return dex.to(self._function_unit, np.log10(x))
/Users/jpierel/CodeBase/jhat/jhat/simple_jwst_phot.py:897: RuntimeWarning: invalid value encountered in log10
phot['magerr'] = 2.5 * np.log10(1.0 + (fluxerr/flux))
aper_sum_4.1px annulus_median_4.1px aper_bkg_4.1px ... x_idl y_idl
-------------- -------------------- -------------- ... ---------- ----------
297.650667 -99.99 -5238.189153 ... 37.152004 -55.873059
378.648855 5.548146 290.651448 ... 30.043803 -54.923731
115.583682 -99.99 -5238.189153 ... -76.001952 -54.071845
294.794123 0.126702 6.63754 ... -76.013673 -51.700329
234.943054 0.218582 11.450907 ... -76.018672 -50.556787
467.419303 6.36806 333.604394 ... -3.833359 -48.74132
342.922546 0.173903 9.110269 ... -76.032743 -47.064001
371.015027 0.207516 10.871191 ... -76.037395 -45.257877
173.426586 0.099091 5.191109 ... -76.03994 -44.794514
838.584206 7.751832 406.096233 ... -9.072841 -43.561009
... ... ... ... ... ...
133.430719 0.1092 5.720696 ... -75.754248 30.346453
245.152683 0.405534 21.244786 ... -75.721452 34.663273
233.690538 0.556301 29.143019 ... -75.716036 35.792661
445.263408 0.72363 37.908921 ... -75.646588 46.09179
176.538 0.783822 41.062174 ... -75.634492 47.785745
617.329614 4.899353 256.663045 ... 36.555135 50.32043
622.234759 4.964599 260.081101 ... 36.55006 51.069329
330.330828 0.707696 37.074164 ... -75.607653 52.275188
378.108748 0.679931 35.619658 ... -75.600128 53.678531
547.845954 -99.99 -5238.189153 ... -75.589207 55.779377
459.438067 -99.99 -5238.189153 ... 36.520526 55.674128
Length = 211 rows
Align the second image
The plots outputted here show the various steps used by jhat to determine the true matching sources in the image, and the subsequent correction needed for optimal alignment.
wcs_align = st_wcs_align()
wcs_align.run_all(align_image,
telescope='jwst',
outsubdir='mastDownload',
refcat_racol='ra',
refcat_deccol='dec',
refcat_magcol='mag',
refcat_magerrcol='dmag',
overwrite=True,
d2d_max=1,
showplots=2,
refcatname=ref_catname,
histocut_order='dxdy',
sharpness_lim=(0.3,0.9),
roundness1_lim=(-0.7, 0.7),
SNR_min= 3,
dmag_max=1.0,
objmag_lim =(14,24))
.. rst-class:: sphx-glr-horizontal
*
.. image-sg:: /examples/images/sphx_glr_plot_c_miri_003.png
:alt: Initial cut: d2d_max=1, dmag_max=None, Nbright=None, delta_mag_lim=(None, None)
:srcset: /examples/images/sphx_glr_plot_c_miri_003.png
:class: sphx-glr-multi-img
*
.. image-sg:: /examples/images/sphx_glr_plot_c_miri_004.png
:alt: dx, dx, dx, slope:4.882812499998248e-05, 3-sigma cut: 72 out of 75 left mean = -0.158 px, stdev = 0.191 px
:srcset: /examples/images/sphx_glr_plot_c_miri_004.png
:class: sphx-glr-multi-img
*
.. image-sg:: /examples/images/sphx_glr_plot_c_miri_005.png
:alt: dy, dy, dy, slope:4.882812499998248e-05, 3-sigma cut: 61 out of 61 left mean = 0.127 px, stdev = 0.257 px
:srcset: /examples/images/sphx_glr_plot_c_miri_005.png
:class: sphx-glr-multi-img
.. rst-class:: sphx-glr-script-out
.. code-block:: none
0 ./mastDownload/jw02107038001_02105_00001_mirimage.phot.txt
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/astropy/wcs/wcs.py:725: FITSFixedWarning: 'datfix' made the change 'Set DATE-BEG to '2022-07-06T17:47:53.158' from MJD-BEG.
Set DATE-AVG to '2022-07-06T17:48:32.008' from MJD-AVG.
Set DATE-END to '2022-07-06T17:49:10.859' from MJD-END'.
warnings.warn(
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/astropy/wcs/wcs.py:725: FITSFixedWarning: 'obsfix' made the change 'Set OBSGEO-L to -72.174733 from OBSGEO-[XYZ].
Set OBSGEO-B to -38.353284 from OBSGEO-[XYZ].
Set OBSGEO-H to 1740817774.322 from OBSGEO-[XYZ]'.
warnings.warn(
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/astropy/stats/sigma_clipping.py:411: AstropyUserWarning: Input data contains invalid values (NaNs or infs), which were automatically clipped.
warnings.warn('Input data contains invalid values (NaNs or '
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/astropy/stats/sigma_clipping.py:411: AstropyUserWarning: Input data contains invalid values (NaNs or infs), which were automatically clipped.
warnings.warn('Input data contains invalid values (NaNs or '
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/astropy/units/function/logarithmic.py:47: RuntimeWarning: invalid value encountered in log10
return dex.to(self._function_unit, np.log10(x))
/Users/jpierel/CodeBase/jhat/jhat/simple_jwst_phot.py:897: RuntimeWarning: invalid value encountered in log10
phot['magerr'] = 2.5 * np.log10(1.0 + (fluxerr/flux))
*** Note: close plot to continue!
slope intercept maxval index d_bestguess fwhm multimax
0.000049 -0.025 54.211365 28 -0.241148 1.0 False
Keeping 75 out of 75, skippin 0 because of null values in columns d_rot_tmp
median: -0.161845
75.000000 percentile cut: max residual for cut: 0.212953
median: -0.162629
i:00 mean:-0.162629(0.014451) stdev:0.107171(0.010127) X2norm:0.99 Nchanged:0 Ngood:56 Nclip:19
mean: -0.144962
i:01 mean:-0.144962(0.017233) stdev:0.136786(0.012090) X2norm:1.00 Nchanged:8 Ngood:64 Nclip:11
mean: -0.145593
i:02 mean:-0.145593(0.019481) stdev:0.159460(0.013674) X2norm:1.00 Nchanged:4 Ngood:68 Nclip:7
mean: -0.158406
i:03 mean:-0.158406(0.020980) stdev:0.174271(0.014729) X2norm:1.00 Nchanged:2 Ngood:70 Nclip:5
mean: -0.158419
i:04 mean:-0.158419(0.022668) stdev:0.191002(0.015917) X2norm:1.00 Nchanged:2 Ngood:72 Nclip:3
mean: -0.158419
i:05 mean:-0.158419(0.022668) stdev:0.191002(0.015917) X2norm:1.00 Nchanged:0 Ngood:72 Nclip:3
slope intercept maxval index d_bestguess fwhm multimax
0.000049 -0.025195 41.262524 36 0.051438 1.0 False
Keeping 61 out of 61, skippin 0 because of null values in columns d_rot_tmp
median: 0.109562
75.000000 percentile cut: max residual for cut: 0.288823
median: 0.102123
i:00 mean:0.102123(0.021089) stdev:0.139885(0.014745) X2norm:0.99 Nchanged:0 Ngood:45 Nclip:16
mean: 0.114686
i:01 mean:0.114686(0.024837) stdev:0.177373(0.017393) X2norm:1.00 Nchanged:7 Ngood:52 Nclip:9
mean: 0.105753
i:02 mean:0.105753(0.028650) stdev:0.214397(0.020080) X2norm:1.00 Nchanged:5 Ngood:57 Nclip:4
mean: 0.126734
i:03 mean:0.126734(0.033214) stdev:0.257275(0.023293) X2norm:1.00 Nchanged:4 Ngood:61 Nclip:0
mean: 0.126734
i:04 mean:0.126734(0.033214) stdev:0.257275(0.023293) X2norm:1.00 Nchanged:0 Ngood:61 Nclip:0
*** Note: close plots to continue!
/Users/jpierel/CodeBase/tweakreg_hack/tweakreg_hack/tweakreg_step_hack.py:540: AstropyDeprecationWarning: The JWSTgWCS class is deprecated and may be removed in a future version.
Use JWSTWCSCorrector instead.
im = JWSTgWCS(
replacing SIP ./mastDownload/jw02107038001_02105_00001_mirimage_jhat.fits
./mastDownload/jw02107038001_02105_00001_mirimage_jhat.fits
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/astropy/wcs/wcs.py:725: FITSFixedWarning: 'obsfix' made the change 'Set OBSGEO-L to -72.174733 from OBSGEO-[XYZ].
Set OBSGEO-B to -38.353284 from OBSGEO-[XYZ].
Set OBSGEO-H to 1740817774.322 from OBSGEO-[XYZ]'.
warnings.warn(
*** Note: close plots to continue!
0
Check the Output
The reference image has not changed, but let's read in the newly aligned image and compare with the original. subsequent correction needed for optimal alignment.
aligned_image = os.path.join('mastDownload',os.path.basename(align_image).replace('cal.fits','jhat.fits'))
aligned_fits = fits.open(aligned_image)
aligned_data = fits.open(aligned_image)['SCI',1].data
aligned_y,aligned_x = skycoord_to_pixel(star_location,wcs.WCS(aligned_fits['SCI',1],aligned_fits))
aligned_cutout = extract_array(aligned_data,(11,11),(aligned_x,aligned_y))
norm3 = simple_norm(aligned_cutout,stretch='log',min_cut=-1,max_cut=200)
fig,axes = plt.subplots(1,3)
axes[0].imshow(ref_cutout, origin='lower',
norm=norm1,cmap='gray')
axes[1].imshow(align_cutout, origin='lower',
norm=norm2,cmap='gray')
axes[2].imshow(aligned_cutout, origin='lower',
norm=norm3,cmap='gray')
axes[0].set_title('Reference')
axes[1].set_title('To Align')
axes[2].set_title('Aligned')
for i in range(3):
axes[i].tick_params(labelcolor='none',axis='both',color='none')
plt.show()
.. image-sg:: /examples/images/sphx_glr_plot_c_miri_006.png :alt: Reference, To Align, Aligned :srcset: /examples/images/sphx_glr_plot_c_miri_006.png :class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
.. code-block:: none
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/astropy/wcs/wcs.py:725: FITSFixedWarning: 'obsfix' made the change 'Set OBSGEO-L to -72.174733 from OBSGEO-[XYZ].
Set OBSGEO-B to -38.353284 from OBSGEO-[XYZ].
Set OBSGEO-H to 1740817774.322 from OBSGEO-[XYZ]'.
warnings.warn(
.. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 32.073 seconds)
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_c_miri.py <plot_c_miri.py>`
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_c_miri.ipynb <plot_c_miri.ipynb>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_