Lecturer: Christoffer Fremling
Jupyter Notebook Authors: Igor Andreoni & Cameron Hummels
This is a Jupyter notebook lesson taken from the GROWTH Winter School 2018. For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-astro-school-2018-resources.html
Learn how to perform image subtraction to discover astronomical transients from multiple consecutive images.
See GROWTH school webpage for detailed instructions on how to install these modules and packages. Nominally, you should be able to install the python modules with pip install <module>
. The external astromatic packages are easiest installed using package managers (e.g., rpm
, apt-get
).
# Import the relevant packages
import numpy as np
from astropy.io import fits #FITS files handling
import os #Call commands from outside Python
from astropy.io import ascii #Read/write ascii files
# Useful to smooth the images with a Gaussian kernel before the subtraction
from scipy.signal import convolve as scipy_convolve
from astropy.convolution import convolve
from astropy.convolution import Gaussian2DKernel
from astropy.stats import sigma_clipped_stats
from astropy.coordinates import SkyCoord
from astropy import units as u
# Plot
import matplotlib.pyplot as plt
# Running external programs
import subprocess
import shutil
In order for this jupyter notebook to function correctly, we must have some external software installed, as described above. The following step assures that these are installed properly before getting to the rest of the content of this lesson.
dependencies = ['SWarp', 'sextractor', 'psfex', 'ds9']
def test_dependency(dep):
try:
subprocess.Popen(dep, stderr=subprocess.PIPE, shell=True).stderr.read()
print("%s is installed properly. OK" % dep)
return 1
except ImportError:
print("===%s IS NOT YET INSTALLED PROPERLY===" % dep)
return 0
i = 0
for dep in dependencies:
i += test_dependency(dep)
print("\n%i out of %i dependencies installed properly." % (i, len(dependencies)))
if i != len(dependencies):
print("Please correctly install these programs before continuing.")
else:
print("You are ready to continue.")
SWarp is installed properly. OK sextractor is installed properly. OK psfex is installed properly. OK ds9 is installed properly. OK 4 out of 4 dependencies installed properly. You are ready to continue.
# remove temporary fits files in current working directory
[os.remove(f) for f in os.listdir() if f.endswith('.fits')]
# Set directory structure
cwd = os.getcwd()
proc_dir = os.path.join(cwd, 'processed')
data_dir = os.path.join(cwd, 'data')
out_dir = os.path.join(proc_dir, 'out')
if os.path.isdir(proc_dir):
shutil.rmtree(proc_dir)
os.mkdir(proc_dir)
for f in os.listdir(data_dir):
shutil.copy2(os.path.join(data_dir, f), os.path.join(proc_dir, f))
os.chdir(proc_dir)
Define the reference and science images. Open them with ds9 to give them a look.
Also, what is the size of the images in pixel? This information will be useful when we want to align the images using Swarp.
# Reference image
ref_image_name = os.path.join(data_dir, 'refimg_i.fits')
ref_image = fits.open(ref_image_name)
#Plot up the reference image
mean, median, std = sigma_clipped_stats(ref_image[0].data)
plt.figure(figsize=(8,8))
#set the scale of the image based on its statistics
plt.imshow(ref_image[0].data, vmin = median - 2*std, vmax = median + 2*std)
plt.colorbar()
plt.title('Reference image')
plt.show()
#Image size?
print("The dimension of the X axis of the reference image is ")
print(ref_image[0].header["NAXIS1"])
print("The dimension of the Y axis of the reference image is ")
print(ref_image[0].header["NAXIS2"])
The dimension of the X axis of the reference image is 1024 The dimension of the Y axis of the reference image is 2048
Let's do the same for the science image. Can you already spot the Supernova?
#Science image
sci_image_name=os.path.join(data_dir, '20120419094409p.fits')
sci_image = fits.open(sci_image_name)
#Open the images
##os.system('ds9 -zscale '+sci_image_name +' ' + ref_image_name +'&')
#Plot up the science image
mean, median, std = sigma_clipped_stats(sci_image[0].data)
plt.figure(figsize=(8,8))
#set the scale of the image based on its statistics
plt.imshow(sci_image[0].data, vmin = median - 2*std, vmax = median + 2*std)
plt.colorbar()
plt.title('Science image')
plt.show()
#Image size?
print("The dimension of the X axis of the science image is ")
print(sci_image[0].header["NAXIS1"])
print("The dimension of the Y axis of the science image is ")
print(sci_image[0].header["NAXIS2"])
The dimension of the X axis of the science image is 1024 The dimension of the Y axis of the science image is 2048
Use the AstrOmatic Swarp package to align the images. Swarp relies on the astrometric information of the image (in other words, on the sky coordinates), therefore both the science and reference images must be astrometrically calibrated (for example, using the AstrOmatic SCAMP package). In this module we assume that the input images are already calibrated.
#Swarp command
try:
command = "SWarp %s %s -c %s -SUBTRACT_BACK N -RESAMPLE Y -RESAMPLE_DIR . -COMBINE N -IMAGE_SIZE 1800,900" % (sci_image_name, ref_image_name, os.path.join(data_dir, 'config.swarp'))
print('Executing command: %s' % command)
rval = subprocess.run(command.split(), check=True)
print('Success!')
except subprocess.CalledProcessError as err:
print('Could not run SWarp with exit error %s'%err)
#Names of the aligned images
sci_image_aligned_name=sci_image_name.replace(".fits", ".resamp.fits").replace('data','processed')
ref_image_aligned_name=ref_image_name.replace(".fits", ".resamp.fits").replace('data','processed')
Executing command: SWarp /home/chummels/new_version/image_subtraction/data/20120419094409p.fits /home/chummels/new_version/image_subtraction/data/refimg_i.fits -c /home/chummels/new_version/image_subtraction/data/config.swarp -SUBTRACT_BACK N -RESAMPLE Y -RESAMPLE_DIR . -COMBINE N -IMAGE_SIZE 1800,900 Success!
If we attempt an image subtraction now, what does the result look like?
#Test image subtraction:
ref_image_aligned=fits.open(ref_image_aligned_name)
hdr_ref=ref_image_aligned[0].header #save fits header
sci_image_aligned=fits.open(sci_image_aligned_name)
hdr_sci=sci_image_aligned[0].header #save fits header
image_sub = sci_image_aligned[0].data-ref_image_aligned[0].data
hdu_image_sub=fits.PrimaryHDU(image_sub)
hdu_image_sub.writeto("sub_test_0.fits")
#Plot up the result of the image subtraction
mean, median, std = sigma_clipped_stats(hdu_image_sub.data)
plt.figure(figsize=(8,8))
#set the scale of the image based on its statistics
plt.imshow(hdu_image_sub.data, vmin = median - 2*std, vmax = median + 2*std, cmap='gray')
plt.colorbar()
plt.title('Test image subtraction')
plt.show()
# Background subtraction. Import the relevant packages
from astropy.stats import SigmaClip
from photutils import Background2D, MedianBackground
from photutils import make_source_mask
mask_sci = make_source_mask(sci_image_aligned[0].data, snr=2, npixels=3, dilate_size=11)
mask_ref = make_source_mask(ref_image_aligned[0].data, snr=2, npixels=3, dilate_size=11)
sci_image_aligned_name=os.path.join(proc_dir, "bg_sub_sci.fits")
ref_image_aligned_name=os.path.join(proc_dir, "bg_sub_ref.fits")
# remove temporary fits files
if os.path.exists(sci_image_aligned_name): os.remove(sci_image_aligned_name)
if os.path.exists(ref_image_aligned_name): os.remove(ref_image_aligned_name)
sigma_clip = SigmaClip(sigma=3, iters=15) # Sigma clipping
bkg_estimator = MedianBackground()
bkg_sci = Background2D(sci_image_aligned[0].data, (200, 150), filter_size=(3, 3), sigma_clip=sigma_clip, bkg_estimator=bkg_estimator, mask=mask_sci)
bkg_ref = Background2D(ref_image_aligned[0].data, (200, 150), filter_size=(3, 3), sigma_clip=sigma_clip, bkg_estimator=bkg_estimator, mask=mask_ref)
#Remove the background from the science image
sci_image_aligned[0].data = sci_image_aligned[0].data-bkg_sci.background
hdu_image_sub=fits.PrimaryHDU(sci_image_aligned[0].data-bkg_sci.background)
hdu_image_sub.writeto(sci_image_aligned_name)
#Remove the background from the reference image
ref_image_aligned[0].data = ref_image_aligned[0].data-bkg_ref.background
hdu_image_sub=fits.PrimaryHDU(ref_image_aligned[0].data-bkg_ref.background)
hdu_image_sub.writeto(ref_image_aligned_name)
WARNING: AstropyDeprecationWarning: "iters" was deprecated in version 3.1 and will be removed in a future version. Use argument "maxiters" instead. [warnings]
What do the background-subtracted images look like?
#Display with ds9
try:
command = "ds9 -zscale %s %s" % (sci_image_aligned_name, ref_image_aligned_name)
print('Executing command: %s' % command)
rval = subprocess.run(command.split(), check=True)
print('Success!')
except subprocess.CalledProcessError as err:
print('Could not run SWarp with exit error %s'%err)
#Plot here the background image
plt.imshow(bkg_sci.background, origin='lower', cmap='Greys_r')
#plt.imshow(sci_image_aligned[0].data-bkg_sci.background, origin='lower', cmap='Greys_r')
Executing command: ds9 -zscale /home/chummels/new_version/image_subtraction/processed/bg_sub_sci.fits /home/chummels/new_version/image_subtraction/processed/bg_sub_ref.fits Success!
<matplotlib.image.AxesImage at 0x7fb7b811e8d0>
The atmosphere heavily affects the PSF of the images by determining the "seeing" conditions. The seeing for ground-based optical telescopes is usually measured as the FWHM of the imaging PSF. Properties of the atmosphere can change very rapidly, so it is rare that science and reference images are characterized by the same seeing. Therefore their PSFs are usually different, which is a problem for image subtraction.
The PSF of the science and reference images can be matched in several different ways. Here we start by performing a first source extraction on both the science image. We can use the catalogs of sources that we obtain for two main purposes:
if os.path.exists('prepsfex.cat'): #Remove possible temporary files
os.remove("prepsfex.cat")
try:
command = "sextractor %s -c %s -CATALOG_NAME %s -MAG_ZEROPOINT 25.0" % (sci_image_aligned_name, os.path.join(data_dir, 'prepsfex.sex'), os.path.join(data_dir, 'prepsfex.cat'))
print('Executing command: %s\n' % command)
rval = subprocess.run(command.split(), check=True)
print('Success!')
except subprocess.CalledProcessError as err:
print('Could not run SExtractor with exit error %s'%err)
Executing command: sextractor /home/chummels/new_version/image_subtraction/processed/bg_sub_sci.fits -c /home/chummels/new_version/image_subtraction/data/prepsfex.sex -CATALOG_NAME /home/chummels/new_version/image_subtraction/data/prepsfex.cat -MAG_ZEROPOINT 25.0 Success!
Now we use another software part of the AstrOmatic suite, PSFex, to measure the PSF of the science image. PSFex estimates the PSF based on the information present in the catalog generated with SExtractor. Then, let's plot the PSF model obtained with PSFex
#Run PSFex to compute PSF, read and display the final model; needs to output to "out" dir.
if not os.path.isdir('out'): os.mkdir('out')
shutil.copy2(os.path.join(data_dir, 'prepsfex.cat'), os.path.join(proc_dir, 'prepsfex.cat'))
os.system("rm proto_prepsfex.fits")
os.system("psfex prepsfex.cat -c psfex_conf.psfex")
psf_sci_image_name=os.path.join(out_dir,'proto_prepsfex.fits')
print(psf_sci_image_name)
psf_sci_image = fits.open(psf_sci_image_name)
plt.imshow(psf_sci_image[0].data[0], cmap='gray')
/home/chummels/new_version/image_subtraction/processed/out/proto_prepsfex.fits
<matplotlib.image.AxesImage at 0x7fb7baba1240>
Now that the kernel is generated, let's convolve the reference image with the PSF of the science frame.
#Convolve the reference image with the PSF of the science frame
if os.path.exists(os.path.join(proc_dir, 'ref_convolved.fits')):
os.remove(os.path.join(proc_dir, 'ref_convolved.fits'))
kernel_sci = psf_sci_image[0].data[0]
ref_image_aligned=fits.open(ref_image_aligned_name)
ref_conv = scipy_convolve(ref_image_aligned[0].data, kernel_sci, mode='same', method='fft')
#Create a new fits file for the convolved image
hdu_ref_conv = fits.PrimaryHDU(ref_conv,hdr_ref)
hdu_ref_conv.writeto(os.path.join(proc_dir, "ref_convolved.fits"))
#Plot up the convolved reference image
mean, median, std = sigma_clipped_stats(hdu_ref_conv.data)
plt.figure(figsize=(8,8))
#set the scale of the image based on its statistics
plt.imshow(hdu_ref_conv.data, vmin = median - 2*std, vmax = median + 2*std)
plt.colorbar()
plt.title('Convolved science image')
plt.show()
Same as above, but this time we generate a kernel with the properties of the PSF of the reference image. Then, we convolve the science image with this kernel.
#SExtractor command for the ref image
if os.path.exists('prepsfex.cat'):
os.remove("prepsfex.cat")
try:
command = command="sextractor %s -c prepsfex.sex -CATALOG_NAME prepsfex.cat -MAG_ZEROPOINT 25.0" % ref_image_aligned_name
print('Executing command: %s\n' % command)
rval = subprocess.run(command.split(), check=True)
print('Success!')
except subprocess.CalledProcessError as err:
print('Could not run SExtractor with exit error %s'%err)
Executing command: sextractor /home/chummels/new_version/image_subtraction/processed/bg_sub_ref.fits -c prepsfex.sex -CATALOG_NAME prepsfex.cat -MAG_ZEROPOINT 25.0 Success!
#Run PSFex to compute PSF, read and display the final model
if os.path.exists(os.path.join(out_dir,'proto_prepsfex.fits')):
os.remove(os.path.join(out_dir, 'proto_prepsfex.fits'))
try:
command = "psfex prepsfex.cat -c psfex_conf.psfex"
print('Executing command: %s\n' % command)
rval = subprocess.run(command.split(), check=True)
print('Success!')
except subprocess.CalledProcessError as err:
print('Could not run SExtractor with exit error %s'%err)
psf_ref_image_name = os.path.join(out_dir, 'proto_prepsfex.fits')
psf_ref_image = fits.open(psf_ref_image_name)
plt.imshow(psf_ref_image[0].data[0], cmap='gray')
Executing command: psfex prepsfex.cat -c psfex_conf.psfex Success!
<matplotlib.image.AxesImage at 0x7fb7ba33a278>
os.system("rm sci_convolved.fits") #Remove possible existing files
kernel_ref = psf_ref_image[0].data[0]
# Read the SCIENCE image and convolve it with the PSF of the reference frame
sci_image_aligned=fits.open(sci_image_aligned_name)
sci_conv = scipy_convolve(sci_image_aligned[0].data, kernel_ref, mode='same', method='fft')
#Create a new fits file for the convolved image
hdu_sci_conv = fits.PrimaryHDU(sci_conv,hdr_sci)
hdu_sci_conv.writeto("sci_convolved.fits")
#Plot up the convolved science image
mean, median, std = sigma_clipped_stats(hdu_sci_conv.data)
plt.figure(figsize=(8,8))
#set the scale of the image based on its statistics
plt.imshow(hdu_sci_conv.data, vmin = median - 2*std, vmax = median + 2*std)
plt.colorbar()
plt.title('Convolved science image')
plt.show()
Now that the science image is convolved with (an approximation of) the PSF of the reference image, and the reference image is convolved with the PSF of the science image, we can perform the image subtraction.
#Test image subtraction:
os.system("rm sub_test.fits")
from image_registration import chi2_shift
from image_registration.fft_tools import shift
import scipy
from scipy import ndimage, misc
import matplotlib.pyplot as plt
import numpy.fft
# Fine tuning of image alignment
xoff, yoff, exoff, eyoff = chi2_shift(ref_image_aligned[0].data, sci_conv, 10, return_error=True, upsample_factor='auto')
print(xoff,yoff)
sci_conv_shift=scipy.ndimage.shift(sci_conv, [-yoff, -xoff], order=3, mode='reflect', cval=0.0, prefilter=True)
/home/chummels/miniconda3/lib/python3.7/site-packages/image_registration/fft_tools/convolve_nd.py:269: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result. bigarray[arrayslices] = array /home/chummels/miniconda3/lib/python3.7/site-packages/image_registration/fft_tools/convolve_nd.py:270: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result. bigkernel[kernslices] = kernel /home/chummels/miniconda3/lib/python3.7/site-packages/image_registration/fft_tools/convolve_nd.py:325: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result. result = rifft[arrayslices].real /home/chummels/miniconda3/lib/python3.7/site-packages/image_registration/fft_tools/zoom.py:101: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result. outarr[ii] = outarr_d[dims]
0.572265625 0.076171875
The science and reference images are usually obtained with different exposure times. In addition, the reference image can be the stack of several images to increase the depth. Finally, different CCDs of the same camera (or even different regions of the same CCD when multiple amplifiers are present) may have slightly different gain.
The background subtraction should have removed the non-linear offsets between science and reference images. We can therefore normalize the two images by computing the ratio of bright star fluxes in the two images. Once again, we use SExtractor to extract the flux and other quantities.
# Normalize images using the stars in the image
os.system("rm prepsfex.cat")
sextractor_command="sextractor sci_convolved.fits -c prepsfex.sex -CATALOG_NAME sci_match.cat -MAG_ZEROPOINT 25.0 -CATALOG_TYPE=ASCII_HEAD"
#Run SExtractor on the science image
os.system(sextractor_command)
cat_sci=ascii.read('sci_match.cat')
os.system("rm prepsfex.cat") #rm temporary file
sextractor_command="sextractor ref_convolved.fits -c prepsfex.sex -CATALOG_NAME ref_match.cat -MAG_ZEROPOINT 25.0 -CATALOG_TYPE=ASCII_HEAD"
#Run SExtractor on the reference image
os.system(sextractor_command)
cat_ref=ascii.read('ref_match.cat')
Match the catalog of sources of the reference and science images. Calculate the ratio between the flux of source in the science image over the flux of sources in the reference image.
c1 = SkyCoord(ra=cat_sci['X_WORLD'], dec=cat_sci['Y_WORLD'])
c2 = SkyCoord(ra=cat_ref['X_WORLD'], dec=cat_ref['Y_WORLD'])
#print(c1)
idx, d2d, d3d = c1.match_to_catalog_3d(c2)
index_arr=[]
ratio_arr=[]
for i, i2, d in zip(idx, np.arange(len(d2d)),d2d):
#print(i,d)
index_arr.append(i)
print(cat_ref['X_IMAGE'][i],cat_ref['Y_IMAGE'][i],' ', cat_sci['X_IMAGE'][i2],cat_sci['Y_IMAGE'][i2])
print(cat_ref['FLUX_AUTO'][i], cat_sci['FLUX_AUTO'][i2] )
ratio_arr.append(cat_sci['FLUX_AUTO'][i2] / cat_ref['FLUX_AUTO'][i])
1307.4404 336.0891 1308.3041 334.046 651472.8 2253524.0 1715.7888 50.3508 1716.6487 50.7292 2919.975 9292.723 989.4865 36.1106 989.7034 36.4439 496.8632 1760.206 824.3469 6.8587 825.5301 6.638 111479.9 320043.0 1326.342 20.3914 1327.2452 20.7808 4020.154 12751.94 486.3967 21.6124 487.562 20.3622 587.6872 2021.381 1004.6854 868.2502 1004.7606 869.0113 627.3826 2062.861 1369.2184 895.2587 1369.7913 893.0807 26483.88 31209.88 1139.5363 818.4732 1153.4922 873.0015 13604.51 627.0291 1534.5369 857.9026 1535.7301 859.6194 728.8522 3655.551 101.7102 857.9676 100.7837 858.4355 1282.349 4143.87 590.6655 814.0704 590.8465 813.879 156837.7 528038.0 704.8513 855.9346 704.2143 855.5515 1331.656 4858.137 942.095 838.1228 941.8206 837.4998 569.4036 1835.311 1094.1759 852.9767 1094.6216 853.0528 3291.538 11230.72 1139.5363 818.4732 1139.7107 818.1718 13604.51 47744.27 26.1865 823.9939 26.0088 824.3149 963.9858 3500.531 1096.5757 820.5004 1095.8405 819.9836 549.3563 3515.232 865.2557 801.5339 865.4846 802.5891 2189.82 7305.347 348.6858 800.5801 348.6887 800.9792 1743.578 6068.929 967.2813 788.6882 966.2339 788.3749 1957.195 6165.656 1526.4729 789.1428 1527.2566 788.9078 1978.886 6543.973 58.4631 792.6835 57.9063 793.7729 339.4756 1844.107 426.8611 785.0707 428.2927 786.1447 958.0468 3790.296 493.6175 761.9312 471.8565 791.9683 580.2625 986.5802 843.3252 771.2743 843.5197 770.931 3493.601 11418.51 70.2446 759.2386 69.329 758.8787 1542.78 6586.326 493.6175 761.9312 493.6544 761.6369 580.2625 1749.379 935.9514 752.2814 935.8901 751.705 2032.108 8306.979 1731.924 752.2356 1732.6384 754.4844 403.6363 952.634 1355.9746 748.3876 1357.3074 748.3896 441.4305 1251.326 1162.9747 745.5599 1164.0452 745.1471 706.9064 3135.146 352.8905 742.475 351.8022 742.4783 373.2087 1922.609 230.4661 739.9128 230.9833 739.1942 536.1561 1118.231 1668.2959 726.557 1668.6002 726.5767 5873.38 21015.47 255.4937 730.6091 254.9797 730.6711 1032.282 2525.637 1794.6632 729.9658 1793.9891 731.6024 272.3101 1214.481 931.6357 725.4043 932.2894 725.1449 557.0006 2520.398 449.1419 716.8032 449.5049 717.9985 165.1539 747.1612 437.221 692.4268 438.67 692.7363 305.4562 1503.349 723.5328 626.7635 723.7566 626.8039 79622.31 274897.4 462.6985 651.155 461.6394 651.0398 244.674 545.7125 278.7191 636.4144 278.6997 636.6849 1314.119 3250.896 211.1954 628.006 211.6926 627.8864 1058.741 3919.839 129.7976 624.9365 129.8166 625.4357 2044.77 7068.42 1342.7594 606.9783 1343.1217 607.4101 275.7008 1711.297 1089.3788 602.7543 1089.5057 602.256 3067.968 10063.96 129.7976 624.9365 103.9308 606.9385 2044.77 1093.857 568.5526 580.5348 568.993 580.5886 1752.51 6159.887 411.4061 574.5271 411.7411 573.8314 989.4938 2389.868 1081.5079 560.5482 1082.2012 560.3478 14970.03 51515.83 1429.6876 577.2727 1430.306 577.1677 3327.188 11034.56 707.9882 890.5663 708.514 890.0958 506.3508 1630.051 855.0543 507.9171 856.2205 506.8423 411.5435 1009.149 278.3008 429.6241 195.9658 494.5762 16498.26 1038.307 350.217 456.4214 350.7393 456.2532 2544.63 9481.262 335.6344 454.9604 335.9444 454.575 2422.551 9658.513 1574.235 474.8135 1575.0111 474.7204 1446.555 5143.764 1717.3147 467.6023 1718.5717 467.4328 5636.654 19934.27 374.2676 541.7355 374.2659 541.265 353.256 613.7535 278.3008 429.6241 278.6964 429.3744 16498.26 57589.93 1463.5232 511.6817 1463.9884 511.82 2818.144 10334.56 433.5878 397.8126 434.0927 397.7597 3460.077 12372.07 552.1423 399.7553 552.3612 398.7852 312.2324 1479.695 644.9651 389.4753 646.224 389.5564 791.691 2884.686 22.5811 388.9146 22.7882 388.811 797.3258 2668.008 1767.859 379.3983 1768.2325 379.8269 1217.849 3842.359 22.5811 388.9146 77.4794 384.1741 797.3258 882.8976 708.9019 379.0284 709.8181 378.7561 348.3668 1664.337 609.8321 373.9482 610.4021 372.6742 703.8999 2877.248 494.4212 358.6339 495.2029 358.3848 5629.839 19451.8 865.9969 359.9691 868.3767 361.2997 496.1679 2247.335 438.2054 345.4932 440.0637 346.7479 333.8689 1189.671 49.1226 314.4248 49.0457 314.6785 427.1208 1821.7 1636.6123 309.1001 1637.4702 309.7581 300.389 784.697 214.6552 282.3522 215.2468 281.8102 80881.34 275954.8 534.7634 296.5371 535.7922 295.7406 465.4218 2215.764 1145.457 224.9014 1070.3884 278.3593 267.4777 1003.717 130.6825 275.8024 131.3335 275.666 442.1817 2657.482 86.5915 268.3777 85.6523 269.4825 485.4805 2446.421 156.5927 265.2281 158.6548 262.7082 271.3367 7655.383 628.6154 228.5701 684.1796 263.2155 375.4998 1756.098 348.2822 257.4895 349.3839 256.7439 400.1925 6949.453 348.2822 257.4895 409.4536 252.0032 400.1925 4581.044 861.645 233.2094 904.9395 266.8738 1313.643 1898.477 348.2822 257.4895 431.0288 244.0005 400.1925 2283.784 1639.571 235.3087 1639.4666 235.1931 703.1941 2224.841 841.9728 228.2311 842.8556 227.6655 3862.442 12282.8 861.645 233.2094 863.1032 232.3071 1313.643 2661.551 100.7696 250.7178 101.9812 250.2885 1143.564 8902.272 1145.457 224.9014 1145.8398 225.355 267.4777 1094.462 1183.364 222.9198 1185.0245 222.8451 223.1765 1137.391 528.9729 219.6311 530.4017 215.5882 904.0701 1127.782 788.6674 196.7207 789.3939 196.3792 620.7944 2055.756 788.6674 196.7207 741.4003 197.6038 620.7944 452.1226 853.0955 190.7964 854.1906 192.1039 668.8994 1159.535 631.7028 190.8574 632.0488 191.1087 307.3482 1544.053 788.6674 196.7207 818.1182 192.19 620.7944 551.9327 198.8466 180.0815 199.8305 180.0823 826.8592 2737.341 24.7322 183.3822 23.7821 181.6046 1394.893 1799.682 323.8443 167.8755 325.0917 167.3712 540.4866 1758.589 1053.9828 134.1353 1055.1072 134.462 1384.323 3923.667 120.0396 131.1292 120.2166 130.3288 284.0773 1197.733 1670.4735 117.4928 1671.1627 117.587 436.7834 1038.625 716.8115 116.4918 717.2381 115.7344 323.8334 1075.751 988.972 104.5019 989.9993 105.2725 571.8413 1518.872 105.7238 81.562 106.636 81.4607 810.347 3038.781 1132.5763 74.0582 1134.6237 73.889 270.9677 825.7562 323.8443 167.8755 341.3011 61.7882 540.4866 1441.423 1513.7533 878.9217 1513.5853 879.0612 4269.81 14705.03 105.7238 81.562 45.8848 40.7333 810.347 1109.834 348.6858 800.5801 242.7746 874.41 1743.578 1798.471 1054.734 866.157 1054.8322 866.0277 2044.373 6996.284 716.8115 116.4918 632.5135 36.3446 323.8334 833.5007 1776.1282 854.8127 1776.3606 854.9861 8772.116 28134.03 177.7335 104.5981 206.5004 34.9997 224.7884 278.9173
Find the scaling factor
scale=np.median(ratio_arr)
print("The scale is", scale)
The scale is 3.417106350177123
Rescale the science image and perform the image subtraction.
os.system("rm sub_test.fits")
# Perform subtraction
image_sub = sci_conv_shift-ref_conv*scale
hdu_image_sub=fits.PrimaryHDU(image_sub)
hdu_image_sub.writeto("sub_final.fits")
Plot the results...
#Plot up the result of the image subtraction
mean, median, std = sigma_clipped_stats(hdu_image_sub.data)
plt.figure(figsize=(8,8))
#set the scale of the image based on its statistics
plt.imshow(hdu_image_sub.data[200:600,1100:1600], vmin = median - 2*std, vmax = median + 2*std, cmap='gray')
plt.plot([280,310],[115,115], "r-" )
plt.plot([258,258],[95,65], "r-" )
plt.colorbar()
plt.title('Final image subtraction')
plt.show()
#...and plot up the same region of sky of science and template images
mean, median, std = sigma_clipped_stats(sci_conv_shift)
plt.figure(figsize=(8,8))
#set the scale of the image based on its statistics
plt.imshow(sci_conv_shift[200:600,1100:1600], vmin = median - 3*std, vmax = median + 3*std, cmap='gray')
plt.plot([280,310],[115,115], "r-" )
plt.plot([258,258],[95,65], "r-" )
plt.colorbar()
plt.title('Science image')
plt.show()
mean, median, std = sigma_clipped_stats(ref_conv)
plt.figure(figsize=(8,8))
#set the scale of the image based on its statistics
plt.imshow(ref_conv[200:600,1100:1600], vmin = median - 3*std, vmax = median + 3*std, cmap='gray')
plt.plot([280,310],[115,115], "r-" )
plt.plot([258,258],[95,65], "r-" )
plt.colorbar()
plt.title('Reference image')
plt.show()
#ds9 visualization
#os.system('ds9 -zscale '+sci_image_aligned_name +' ' + ref_image_aligned_name + ' sub_final.fits &')
There you go, now you should be able to easily spot the supernova as a white excess on the grey background!
Note that the bright center of the host galaxy was not perfectly subtracted and left a spurious signal that could be mistaken for real luminosity variability.
These operations can be made automatic and can be incorporated in pipelines that discover transients. Moreover, using the methods learnt in the Photometry module, you can perform forced PSF photometry on the image subtraction to obtain flux measurement of the transient free from the host galaxy contamination.
There are some packages that perform most of the operations above automatically. One of the most popular is called HOTPANTS, which stands for "High Order Transform of PSF ANd Template Subtraction".
The code can be found on GitHub: https://github.com/acbecker/hotpants
To proceed, install this program and continue.
#hotpants
os.system("path/to/hotpants -inim"+" "+sci_image_aligned_name+" "+"-tmplim "+ref_image_aligned_name+" -outim out.fits -oci oci.fits -tg 1 -tr 1 -ig 1 -ir 1 -nrx 1 -nry 1 -nsx 2 -nsy 2 -ng 3 6 0.70 4 1.50 2 3.00 -v 0 -bgo 0 -tl 1 -tu 150000 -il 1 -iu 150000 -ko 0")
hotpants_image = fits.open('oci.fits')
os.system("rm sub_test.fits")
image_sub = hotpants_image[0].data-ref_image_aligned[0].data*0.975
hdu_image_sub=fits.PrimaryHDU(image_sub)
hdu_image_sub.writeto("sub_test_hp.fits")
os.system('ds9 -zscale sub_test_hp.fits &')
We can also perform image subtraction after degrading (i.e., convolving with a kernel) only one of the two images. In fact, in some cases it is preferred to keep the science image untouched and manipulate only the reference to avoid introducing errors in photometry measurements. On the other hand, convolving both images usually leads to "cleaner" image subtraction.
Below are a list of commands to do a better job with manipulating only one image usinf fft techniques. In this example, we manipulate the science image instead of the reference simply because of what the specific characteristics of the images that we are playing with. Identical methods can be applied to manipulating reference images only.
# alternatively to degrade one image only
Im1_psf=np.fft.fftshift(kernel_ref);
Im2_psf=np.fft.fftshift(kernel_sci);
psf_fft1=np.fft.fft2(Im1_psf);
psf_fft2=np.fft.fft2(Im2_psf);
kern=psf_fft1/psf_fft2;
kern_m=np.real(np.fft.fftshift(np.fft.ifft2(kern)))
kern_mean=np.mean(kern_m)
psf_m=(kern_m/(kern_mean*59*59))
plt.imshow(psf_m, cmap='gray')
sci_conv = scipy_convolve(sci_image_aligned[0].data, psf_m, mode='same', method='fft')
from image_registration import chi2_shift
from image_registration.fft_tools import shift
import scipy
from scipy import ndimage, misc
import matplotlib.pyplot as plt
import numpy.fft
xoff, yoff, exoff, eyoff = chi2_shift(ref_image_aligned[0].data, sci_conv, 10, return_error=True, upsample_factor='auto')
print(xoff,yoff)
sci_conv_shift=scipy.ndimage.shift(sci_conv, [-yoff, -xoff], order=3, mode='reflect', cval=0.0, prefilter=True)
os.system("rm sub_test2.fits")
image_sub = sci_conv_shift*scale-ref_image_aligned[0].data
hdu_image_sub=fits.PrimaryHDU(image_sub)
hdu_image_sub.writeto("sub_test2.fits")
os.system('ds9 -zscale sub_test2.fits&')