Image Subtraction Module

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

Objective

Learn how to perform image subtraction to discover astronomical transients from multiple consecutive images.

Key steps

  • Register science and reference images (make them aligned and of the same size)
  • PSF extraction, using PSFex
  • PSF matching by convolution
  • [Zero-point calibration]
  • Image subtraction

Required dependencies

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).

Python modules

  • python 3
  • astropy
  • numpy
  • scipy
  • matplotlib
  • pytest

External packages

In [1]:
# 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

Test dependencies

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.

In [2]:
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.

Set up paths and clear old temp files

In [3]:
# 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)

Reference and science images

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.

In [4]:
# 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?

In [5]:
#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

Align the images

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.

In [6]:
#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?

In [7]:
#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

  • Mask sources in images
  • Use 3 sigma clipping, 15 iterations to filter data and accurately measure the backgorund
  • Then split image into 300x300 pixel boxes and apply 2x2 median filter
In [8]:
# 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?

In [9]:
#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!
Out[9]:
<matplotlib.image.AxesImage at 0x7fb7b811e8d0>