Note
Click here to download the full example code
JWST NIRCAM
Aligning JWST/NIRCAM images with JHAT.
An example NIRCam Dataset is downloaded, and then a series of alignment methods are used. For more information on the key parameters used for alignment see 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
Relative Alignment
Download some Data
For this example we download 2 JWST NIRCam images from MAST. They’re the same field but different filters. Note that the code will also work for level 3 data images.
obs_table1 = Observations.query_criteria(obs_id='jw02107-o041_t019_nircam_clear-f200w')
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-o041_t019_nircam_clear-f360m')
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')
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:JWST/product/jw02107041001_02101_00001_nrcb1_cal.fits to ./mastDownload/JWST/jw02107041001_02101_00001_nrcb1/jw02107041001_02101_00001_nrcb1_cal.fits ... [Done]
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:JWST/product/jw02107041001_02101_00001_nrcblong_cal.fits to ./mastDownload/JWST/jw02107041001_02101_00001_nrcblong/jw02107041001_02101_00001_nrcblong_cal.fits ... [Done]
Examine the Reference Image
ref_image = glob.glob('mastDownload/JWST/*nrcb1*/*cal.fits')[0]
ref_fits = fits.open(ref_image)
ref_data = fits.open(ref_image)['SCI',1].data
norm1 = simple_norm(ref_data,stretch='linear',min_cut=-.5,max_cut=3)
plt.imshow(ref_data, origin='lower',
norm=norm1,cmap='gray')
plt.gca().tick_params(labelcolor='none',axis='both',color='none')
plt.show()
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:41.0532','-43:26:41.128',unit=(u.hourangle,u.deg))
align_image = glob.glob('mastDownload/JWST/*long*/*cal.fits')[0]
align_fits = fits.open(align_image)
align_fits['SCI',1].header['CRPIX1']+=1
align_fits['SCI',1].header['CRPIX2']+=1
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='linear',min_cut=-.5,max_cut=3)
norm2 = simple_norm(align_cutout,stretch='linear',min_cut=-.5,max_cut=3)
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()
/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-06T19:16:42.721' from MJD-BEG.
Set DATE-AVG to '2022-07-06T19:17:14.932' from MJD-AVG.
Set DATE-END to '2022-07-06T19:17:47.142' 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.164999 from OBSGEO-[XYZ].
Set OBSGEO-B to -38.353872 from OBSGEO-[XYZ].
Set OBSGEO-H to 1740894174.999 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-06T19:16:42.721' from MJD-BEG.
Set DATE-AVG to '2022-07-06T19:17:14.932' from MJD-AVG.
Set DATE-END to '2022-07-06T19:17:47.142' 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.164999 from OBSGEO-[XYZ].
Set OBSGEO-B to -38.353872 from OBSGEO-[XYZ].
Set OBSGEO-H to 1740894174.999 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,ee_radius=80)
ref_catname = ref_image.replace('.fits','.phot.txt') # the default
refcat = Table.read(ref_catname,format='ascii')
print(refcat)
0 mastDownload/JWST/jw02107041001_02101_00001_nrcb1/jw02107041001_02101_00001_nrcb1_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-06T19:16:42.721' from MJD-BEG.
Set DATE-AVG to '2022-07-06T19:17:14.932' from MJD-AVG.
Set DATE-END to '2022-07-06T19:17:47.142' 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.164999 from OBSGEO-[XYZ].
Set OBSGEO-B to -38.353872 from OBSGEO-[XYZ].
Set OBSGEO-H to 1740894174.999 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_5.7px annulus_median_5.7px aper_bkg_5.7px ... x_idl y_idl
-------------- -------------------- -------------- ... ---------- ----------
105.310801 0.910453 92.21433 ... -11.710449 -31.373861
121.56116 0.960956 97.329433 ... 1.390701 -31.374766
60.232352 0.498815 50.521961 ... 18.604013 -31.26976
126.737312 1.019729 103.282174 ... -31.055244 -31.32962
113.529229 0.802048 81.234601 ... -18.547597 -31.34987
105.370409 0.767367 77.721958 ... -13.204667 -31.345555
89.411974 0.648643 65.697099 ... 11.290307 -31.300338
110.788413 0.861123 87.217932 ... -0.18106 -31.287703
124.473601 0.862357 87.342937 ... 6.266814 -31.255725
108.728849 0.849587 86.04952 ... 10.634677 -31.259596
... ... ... ... ... ...
33.278942 0.200456 20.302942 ... -18.177622 29.366304
23.6942 0.04802 4.863676 ... -17.540134 29.519926
18.277825 0.013074 1.324192 ... -9.43527 29.52744
0.06818 -99.99 -10127.382506 ... -2.732965 29.758751
14.784146 -99.99 -10127.382506 ... 20.030247 29.877627
47.576291 0.323802 32.795919 ... 8.796703 31.023101
36.125728 0.251559 25.478925 ... 26.612871 31.127561
39.783969 0.184818 18.719089 ... -15.095758 31.1144
341.811179 0.27702 28.0577 ... -9.937289 31.082671
53.600708 0.336015 34.032935 ... 8.547109 31.138805
50.031236 0.275069 27.860125 ... 13.180869 31.166268
Length = 2756 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))
0 ./mastDownload/jw02107041001_02101_00001_nrcblong.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-06T19:16:42.721' from MJD-BEG.
Set DATE-AVG to '2022-07-06T19:17:14.932' from MJD-AVG.
Set DATE-END to '2022-07-06T19:17:47.142' 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.164999 from OBSGEO-[XYZ].
Set OBSGEO-B to -38.353872 from OBSGEO-[XYZ].
Set OBSGEO-H to 1740894174.999 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.000146 -0.15 83.958993 76 -0.50978 0.8 False
Keeping 109 out of 109, skippin 0 because of null values in columns d_rot_tmp
median: -0.572386
75.000000 percentile cut: max residual for cut: 0.188131
median: -0.587588
i:00 mean:-0.587588(0.010388) stdev:0.092914(0.007300) X2norm:0.99 Nchanged:0 Ngood:81 Nclip:28
mean: -0.580012
i:01 mean:-0.580012(0.012575) stdev:0.121920(0.008845) X2norm:1.00 Nchanged:14 Ngood:95 Nclip:14
mean: -0.557890
i:02 mean:-0.557890(0.014267) stdev:0.143378(0.010038) X2norm:1.00 Nchanged:7 Ngood:102 Nclip:7
mean: -0.554064
i:03 mean:-0.554064(0.015446) stdev:0.157514(0.010870) X2norm:1.00 Nchanged:3 Ngood:105 Nclip:4
mean: -0.554064
i:04 mean:-0.554064(0.015446) stdev:0.157514(0.010870) X2norm:1.00 Nchanged:0 Ngood:105 Nclip:4
slope intercept maxval index d_bestguess fwhm multimax
-1.734723e-17 1.776357e-14 79.085569 6 -0.534655 0.8 False
Keeping 103 out of 103, skippin 0 because of null values in columns d_rot_tmp
median: -0.559999
75.000000 percentile cut: max residual for cut: 0.191320
median: -0.559999
i:00 mean:-0.559999(0.010883) stdev:0.094879(0.007646) X2norm:0.99 Nchanged:0 Ngood:77 Nclip:26
mean: -0.554858
i:01 mean:-0.554858(0.013144) stdev:0.124003(0.009243) X2norm:1.00 Nchanged:13 Ngood:90 Nclip:13
mean: -0.565181
i:02 mean:-0.565181(0.015176) stdev:0.148691(0.010675) X2norm:1.00 Nchanged:7 Ngood:97 Nclip:6
mean: -0.564816
i:03 mean:-0.564816(0.016002) stdev:0.158407(0.011258) X2norm:1.00 Nchanged:2 Ngood:99 Nclip:4
mean: -0.569514
i:04 mean:-0.569514(0.016529) stdev:0.164459(0.011629) X2norm:1.00 Nchanged:1 Ngood:100 Nclip:3
mean: -0.569514
i:05 mean:-0.569514(0.016529) stdev:0.164459(0.011629) X2norm:1.00 Nchanged:0 Ngood:100 Nclip:3
*** 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/jw02107041001_02101_00001_nrcblong_jhat.fits
./mastDownload/jw02107041001_02101_00001_nrcblong_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.164999 from OBSGEO-[XYZ].
Set OBSGEO-B to -38.353872 from OBSGEO-[XYZ].
Set OBSGEO-H to 1740894174.999 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='linear',min_cut=-.5,max_cut=3)
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()
/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.164999 from OBSGEO-[XYZ].
Set OBSGEO-B to -38.353872 from OBSGEO-[XYZ].
Set OBSGEO-H to 1740894174.999 from OBSGEO-[XYZ]'.
warnings.warn(
Align to Catalog
You can also align each image to the Gaia DR3 catalog, or you could replace the catalog created in step one with your own catalog of the field.
wcs_align.run_all(align_image,
telescope='jwst',
outsubdir='mastDownload',
overwrite=True,
d2d_max=.5,
showplots=0,
refcatname='Gaia',
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))
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='linear',min_cut=-.5,max_cut=3)
fig,axes = plt.subplots(1,2)
axes[0].imshow(align_cutout, origin='lower',
norm=norm2,cmap='gray')
axes[1].imshow(aligned_cutout, origin='lower',
norm=norm3,cmap='gray')
axes[0].set_title('To Align')
axes[1].set_title('Aligned')
for i in range(2):
axes[i].tick_params(labelcolor='none',axis='both',color='none')
plt.show()
0 ./mastDownload/jw02107041001_02101_00001_nrcblong.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-06T19:16:42.721' from MJD-BEG.
Set DATE-AVG to '2022-07-06T19:17:14.932' from MJD-AVG.
Set DATE-END to '2022-07-06T19:17:47.142' 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.164999 from OBSGEO-[XYZ].
Set OBSGEO-B to -38.353872 from OBSGEO-[XYZ].
Set OBSGEO-H to 1740894174.999 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))
INFO: Query finished. [astroquery.utils.tap.core]
Number of stars: 21
### NO propoer motion correction!!!
/Users/jpierel/miniconda3/envs/tweakreg/lib/python3.10/site-packages/pandas/core/arraylike.py:402: RuntimeWarning: invalid value encountered in sqrt
result = getattr(ufunc, method)(*inputs, **kwargs)
Number of stars after removing nan's: 21
slope intercept maxval index d_bestguess fwhm multimax
0.001367 -1.4 3.661515 35 0.706414 0.8 False
Keeping 5 out of 5, skippin 0 because of null values in columns d_rot_tmp
median: 0.644163
75.000000 percentile cut: max residual for cut: 0.071771
median: 0.646769
i:00 mean:0.646769(0.001950) stdev:0.002758(0.001126) X2norm:0.79 Nchanged:0 Ngood:3 Nclip:2
mean: 0.646657
i:01 mean:0.646657(0.001725) stdev:0.002440(0.000996) X2norm:1.00 Nchanged:0 Ngood:3 Nclip:2
slope intercept maxval index d_bestguess fwhm multimax
0.000537 -0.55 3.0 4 -1.643322 0.8 False
Keeping 3 out of 3, skippin 0 because of null values in columns d_rot_tmp
median: -1.609357
i:00 mean:-1.609357(0.083353) stdev:0.117879(0.048124) X2norm:0.79 Nchanged:0 Ngood:3 Nclip:0
mean: -1.633248
i:01 mean:-1.633248(0.070913) stdev:0.100286(0.040942) X2norm:1.00 Nchanged:0 Ngood:3 Nclip:0
/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/jw02107041001_02101_00001_nrcblong_jhat.fits
./mastDownload/jw02107041001_02101_00001_nrcblong_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.164999 from OBSGEO-[XYZ].
Set OBSGEO-B to -38.353872 from OBSGEO-[XYZ].
Set OBSGEO-H to 1740894174.999 from OBSGEO-[XYZ]'.
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.164999 from OBSGEO-[XYZ].
Set OBSGEO-B to -38.353872 from OBSGEO-[XYZ].
Set OBSGEO-H to 1740894174.999 from OBSGEO-[XYZ]'.
warnings.warn(
Total running time of the script: ( 1 minutes 42.234 seconds)