Tutorial - Part #5 - Image Coaddition

For image subtraction the package has a module called propersubtract, which implements a main diff function.



In [1]:

    
from copy import copy
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors

%matplotlib inline



In [2]:

    
from astropy.visualization import LinearStretch, LogStretch
from astropy.visualization import ZScaleInterval, MinMaxInterval
from astropy.visualization import ImageNormalize



In [3]:

    
palette = copy(plt.cm.gray)
palette.set_bad('r', 1.0)
#palette.set_under('r', 1.0)
#palette.set_over('r', 1.0)



In [4]:

    
import properimage.single_image as si
from properimage import propercoadd as pc



In [5]:

    
im1_path = './../../../data/aligned_eso085-030-004.fit'
im2_path = './../../../data/aligned_eso085-030-005.fit'

To get the subtraction we need to run this function by using both paths for example:



In [6]:

    
from astropy.io.fits import getdata

plt.figure(figsize=(16, 12))
plt.subplot(121)
ref = getdata(im1_path)
norm = ImageNormalize(ref, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.imshow(ref, cmap=plt.cm.gray, norm=norm, interpolation='none')

plt.colorbar(orientation='horizontal')

plt.subplot(122)
ref = getdata(im2_path)
norm = ImageNormalize(ref, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.imshow(ref, cmap=plt.cm.gray, norm=norm, interpolation='none')

plt.colorbar(orientation='horizontal')









    Out[6]:





<matplotlib.colorbar.Colorbar at 0x7f4f957efc90>



In [7]:

    
R, P_r, mask = pc.stack_R(align=False, si_list=[im1_path, im2_path], inf_loss=0.25)









    



Sources found = 84
stamps will be 15 x 15
Sources found = 83
stamps will be 15 x 15
(83, 83, 84)
Matching result::  IDs > 0. => 80
p=2, q=81
('Masked pixels: ', 175)






    



/home/bruno/.virtualenvs/iPTF/local/lib/python2.7/site-packages/astropy/convolution/convolve.py:673: 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/bruno/.virtualenvs/iPTF/local/lib/python2.7/site-packages/astropy/convolution/convolve.py:679: 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] = normalized_kernel
/home/bruno/.virtualenvs/iPTF/local/lib/python2.7/site-packages/astropy/convolution/convolve.py:695: 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.
  bigimwt[arrayslices] = 1.0 - nanmaskarray * interpolate_nan
/home/bruno/.virtualenvs/iPTF/local/lib/python2.7/site-packages/astropy/convolution/convolve.py:703: 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.
  bigimwt[arrayslices] = wtsm.real[arrayslices]
/home/bruno/.virtualenvs/iPTF/local/lib/python2.7/site-packages/astropy/convolution/convolve.py:723: RuntimeWarning: divide by zero encountered in true_divide
  rifft = (ifftn(fftmult)) / bigimwt
/home/bruno/.virtualenvs/iPTF/local/lib/python2.7/site-packages/astropy/convolution/convolve.py:738: 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






    



updating stamp shape to (21,21)
((84, 84), (441, 84))
('Masked pixels: ', 180)
updating stamp shape to (21,21)
((83, 83), (441, 83))
(682, 1024)
starting new process
(682, 1024)
starting new process
all chunks started, and procs appended






    



/home/bruno/.virtualenvs/iPTF/local/lib/python2.7/site-packages/pyfftw/builders/_utils.py:161: 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.
  a_copy[update_input_array_slicer])
/home/bruno/.virtualenvs/iPTF/local/lib/python2.7/site-packages/pyfftw/builders/_utils.py:161: 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.
  a_copy[update_input_array_slicer])






    



loading pickles
loading pickles
S calculated, now starting to join processes
waiting for procs to finish
waiting for procs to finish
processes finished, now returning R

The result is a list of numpy arrays.

The arrays are in order: R, P_r, mask



In [8]:

    
norm = ImageNormalize(R.real, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.figure(figsize=(22, 22))
plt.imshow(np.ma.MaskedArray(R.real, mask=mask), 
           cmap=palette, norm=norm, interpolation='none')

plt.colorbar(orientation='horizontal')









    Out[8]:





<matplotlib.colorbar.Colorbar at 0x7f4f95762910>



In [9]:

    
norm = ImageNormalize(P_r.real, interval=MinMaxInterval(),
                      stretch=LogStretch())
xc, yc = np.where(P_r.real==P_r.real.max())
xc, yc = np.round(xc[0]), np.round(yc[0])

plt.imshow(P_r.real[xc-10:xc+10, yc-10:yc+10], norm=norm,
           cmap='viridis', interpolation='none')
plt.colorbar()









    Out[9]:





<matplotlib.colorbar.Colorbar at 0x7f4f8e40f490>



In [10]:

    
R_1, P_r_1, mask_1 = pc.stack_R(align=False, si_list=[im1_path, im2_path], inf_loss=0.10)









    



Sources found = 84
stamps will be 15 x 15
Sources found = 83
stamps will be 15 x 15
(83, 83, 84)
Matching result::  IDs > 0. => 80
p=2, q=81
('Masked pixels: ', 175)
updating stamp shape to (21,21)
((84, 84), (441, 84))
('Masked pixels: ', 180)
updating stamp shape to (21,21)
((83, 83), (441, 83))
(682, 1024)
starting new process
(682, 1024)
starting new process
all chunks started, and procs appended






    



/home/bruno/.virtualenvs/iPTF/local/lib/python2.7/site-packages/pyfftw/builders/_utils.py:161: 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.
  a_copy[update_input_array_slicer])
/home/bruno/.virtualenvs/iPTF/local/lib/python2.7/site-packages/pyfftw/builders/_utils.py:161: 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.
  a_copy[update_input_array_slicer])






    



loading pickles
loading pickles
S calculated, now starting to join processes
waiting for procs to finish
waiting for procs to finish
processes finished, now returning R

The result is a list of numpy arrays.

The arrays are in order: R, P_r, mask



In [11]:

    
norm = ImageNormalize(R_1.real, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.figure(figsize=(22, 22))
plt.imshow(np.ma.MaskedArray(R_1.real, mask=mask_1), 
           cmap=palette, norm=norm, interpolation='none')

plt.colorbar(orientation='horizontal')









    Out[11]:





<matplotlib.colorbar.Colorbar at 0x7f4f95680190>



In [12]:

    
norm = ImageNormalize(P_r.real, interval=MinMaxInterval(),
                      stretch=LogStretch())
xc, yc = np.where(P_r.real==P_r.real.max())
xc, yc = np.round(xc[0]), np.round(yc[0])
plt.imshow(P_r.real[xc-10:xc+10, yc-10:yc+10], norm=norm,
           cmap='viridis', interpolation='none')
plt.colorbar()









    Out[12]:





<matplotlib.colorbar.Colorbar at 0x7f4f8e2ce4d0>

Coadd using SWarp and compare photometry



In [129]:

    
import shlex, subprocess
cmd = 'swarp -IMAGEOUT_NAME coadd.fits -WEIGHT MAP_VARIANCE '+im1_path+' '+im2_path
cmd = shlex.split(cmd)
subprocess.call(cmd)









    Out[129]:





0



In [130]:

    
plt.figure(figsize=(22, 22))
coadd = getdata('coadd.fits')[18:-18, 28:-28]
norm = ImageNormalize(coadd, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.imshow(coadd, cmap=plt.cm.gray, norm=norm, interpolation='none')

plt.colorbar(orientation='horizontal')









    Out[130]:





<matplotlib.colorbar.Colorbar at 0x7f4f942eff90>

Detect sources on properimage's



In [131]:

    
import sep



In [132]:

    
xc, yc = np.where(P_r.real==P_r.real.max())
xc, yc = np.round(xc[0]), np.round(yc[0])
ker = P_r.real[xc-10:xc+10, yc-10:yc+10]

bkg_p = sep.Background(R.real.astype('<f4'), mask=mask)
srcs_p = sep.extract(R.real.astype('<f4')-bkg_p.back(), 1.5, 
                     err=bkg_p.globalrms, mask=mask)#, filter_kernel=ker)



In [133]:

    
len(srcs_p)









    Out[133]:





861



In [145]:

    
flux_p, fluxerr_p, flag_p = sep.sum_circle(R.real.astype('<f4')-bkg_p.back(), 
                                     srcs_p['x'], srcs_p['y'],
                                     5.0, err=bkg_p.globalrms, gain=1.0, bkgann=(8., 11.))



In [146]:

    
norm = ImageNormalize(R.real, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.figure(figsize=(22, 22))
plt.imshow(np.ma.MaskedArray(R.real, mask=mask), 
           cmap=palette, norm=norm, interpolation='none')
plt.plot(srcs_p['x'], srcs_p['y'], 'bo')
plt.xlim(-20, 1040)
plt.ylim(700, -10)
plt.colorbar(orientation='horizontal')









    Out[146]:





<matplotlib.colorbar.Colorbar at 0x7f4f8d2c5810>



In [147]:

    
bkg_c = sep.Background(coadd.astype('<f4'))
srcs_c = sep.extract(coadd.astype('<f4')-bkg_c.back(), 1.5, err=bkg_c.globalrms)



In [148]:

    
len(srcs_c)









    Out[148]:





949



In [149]:

    
flux, fluxerr, flag = sep.sum_circle(coadd.astype('<f4')-bkg_c.back(), 
                                     srcs_c['x'], srcs_c['y'],
                                     5.0, err=bkg_c.globalrms, gain=1.0, bkgann=(8., 11.))



In [150]:

    
norm = ImageNormalize(coadd, interval=ZScaleInterval(),
                      stretch=LinearStretch())

plt.figure(figsize=(22, 22))
plt.imshow(coadd, cmap=palette, norm=norm, interpolation='none')
plt.plot(srcs_c['x'], srcs_c['y'], 'bo')
plt.xlim(-20, 1040)
plt.ylim(700, -10)
plt.colorbar(orientation='horizontal')









    Out[150]:





<matplotlib.colorbar.Colorbar at 0x7f4f7f9af7d0>



In [157]:

    
mag_p = -2.5*np.log10(flux_p)
mag_c = -2.5*np.log10(flux)









    



/home/bruno/.virtualenvs/iPTF/lib/python2.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in log10
  """Entry point for launching an IPython kernel.
/home/bruno/.virtualenvs/iPTF/lib/python2.7/site-packages/ipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in log10



In [154]:

    
plt.subplot(121)
plt.plot(flux_p, flux_p/fluxerr_p, 'xb')
plt.subplot(122)
plt.plot(flux, flux/fluxerr, 'xr')









    Out[154]:





[<matplotlib.lines.Line2D at 0x7f4f8cb95890>]



In [159]:

    
plt.subplot(121)
plt.plot(mag_p, flux_p/fluxerr_p, 'xb')
plt.subplot(122)
plt.plot(mag_c, flux/fluxerr, 'xr')









    Out[159]:





[<matplotlib.lines.Line2D at 0x7f4f8c90b310>]



In [163]:

    
plt.subplot(121)
plt.hist(mag_p[~np.isnan(mag_p)])
plt.subplot(122)
plt.hist(mag_c[~np.isnan(mag_c)])









    Out[163]:





(array([  7.,  12.,  69., 249., 399., 128.,   8.,   4.,   0.,   1.]),
 array([-16.39415708, -14.67919319, -12.9642293 , -11.24926542,
         -9.53430153,  -7.81933765,  -6.10437376,  -4.38940987,
         -2.67444599,  -0.9594821 ,   0.75548179]),
 <a list of 10 Patch objects>)



In [ ]: