In [8]:
from matplotlib import pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = 12,8
In [9]:
from astropy.io import fits
import photutils
from bok.header_fix import header_fix
from common import *
from bok.constant import const
In [10]:
amp_fits = '/data/d7420.0130.04.fits'
amp_ldac = '/data/d7420.0130.04.ldac'
In [11]:
cata = fits.getdata(amp_ldac, 2)
In [12]:
cata.columns
Out[12]:
In [13]:
c1 = fits.Column(name="MAG_CORR", format="1E")
c2 = fits.Column(name="MAGERR_CORR", format="1E")
cols = fits.ColDefs([ c1, c2 ])
In [63]:
cata.BACKGROUND
Out[63]:
In [17]:
newtab = fits.BinTableHDU.from_columns(cata.columns + cols)
In [20]:
newtab.data.MAG_CORR
Out[20]:
In [14]:
image = fits.getdata(amp_fits) - 32768
In [15]:
plt.imshow(np.log10(image), clim=(0,1))
plt.plot(cata.X_IMAGE, cata.Y_IMAGE, "+")
plt.xlim(10,2000)
plt.ylim(10,2000)
plt.colorbar()
Out[15]:
In [34]:
sum(( (10 < cata.X_IMAGE) & (cata.X_IMAGE < 2000) & (10 < cata.Y_IMAGE) & (cata.Y_IMAGE < 2000) ))
Out[34]:
In [38]:
def noborder(x, y, xlim, ylim) :
return (xlim[0] < x) & (x < xlim[1]) & (ylim[0] < y) & (y < ylim[1])
In [16]:
from astropy.stats import sigma_clipped_stats
mean, median, std = sigma_clipped_stats(image, sigma=3.0, iters=5)
print (mean, median, std)
In [50]:
sigma_clipped_stats(cata.FWHM_IMAGE)
Out[50]:
In [51]:
sdao = photutils.daofind(image, fwhm=4.0, threshold=3.*std)
sdao
Out[51]:
In [52]:
si = photutils.irafstarfind(image, fwhm=4.0, threshold=3.*std)
si
Out[52]:
In [53]:
ixdao = np.where(noborder(sdao["xcentroid"], sdao["ycentroid"], (10, 2000), (10, 2000)))
#ixsou = np.where(noborder(sources["xcentroid"], sources["ycentroid"], (10, 2000), (10, 2000)))
ixi = np.where(noborder(si["xcentroid"], si["ycentroid"], (10, 2000), (10, 2000)))
plt.plot(cata.X_IMAGE, cata.Y_IMAGE, 's')
plt.plot(sdao[ixdao]["xcentroid"], sdao[ixdao]["ycentroid"], '+')
#plt.plot(sources[ixsou]["xcentroid"], sources[ixsou]["ycentroid"], 's')
plt.plot(si[ixi]["xcentroid"], si[ixi]["ycentroid"], 'x')
print (len(ixdao[0]), len(ixi[0]))
In [54]:
from photutils import CircularAperture
In [55]:
pos = list(zip(cata.X_IMAGE, cata.Y_IMAGE))
In [56]:
apertures = CircularAperture(pos, r=3.)
In [58]:
phot_aper_3 = photutils.aperture_photometry(image, apertures)
phot_aper_3
Out[58]:
In [75]:
aper_flux = phot_aper_3["aperture_sum"] - 3*3*np.pi*(cata.BACKGROUND-32768.0)
aper_mag = 25.0 - 2.5 * np.log10(aper_flux)
In [67]:
sigma_clipped_stats(cata.BACKGROUND-32768)
Out[67]:
In [69]:
sigma_clipped_stats(image)
Out[69]:
In [76]:
mag_d = aper_mag - cata.MAG_AUTO
In [83]:
plt.hist(mag_d, range=(-0.5,2.5), bins=30)
Out[83]:
In [ ]: