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]:
In [7]:
R, P_r, mask = pc.stack_R(align=False, si_list=[im1_path, im2_path], inf_loss=0.25)
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]:
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]:
In [10]:
R_1, P_r_1, mask_1 = pc.stack_R(align=False, si_list=[im1_path, im2_path], inf_loss=0.10)
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]:
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]:
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]:
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]:
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]:
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]:
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]:
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]:
In [157]:
mag_p = -2.5*np.log10(flux_p)
mag_c = -2.5*np.log10(flux)
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]:
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]:
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]:
In [ ]: