In [3]:
import os
import sys
import numpy as np
import matplotlib
matplotlib.use('nbagg')
#from matplotlib import style
#style.use('ggplot')
import matplotlib.pyplot as plt
from photutils import Background2D, MedianBackground
import photutils
from scipy import optimize
from astroscrappy import detect_cosmics
from skimage.morphology import reconstruction
from scipy.ndimage.filters import median_filter
from scipy.ndimage import gaussian_filter
from scipy.signal import medfilt
from skimage.filters.rank import median
import astropy.units as u
from astropy.modeling.models import Gaussian2D, Polynomial2D
from astropy.modeling.fitting import LevMarLSQFitter
from astropy.coordinates import SkyCoord, match_coordinates_3d
from astroscrappy import detect_cosmics
%load_ext autoreload
%autoreload 2
from astropy.io import fits
from mmtwfs.wfs import *
from mmtwfs.zernike import ZernikeVector
In [4]:
def implot(im):
norm = wfs_norm(im)
plt.imshow(im, norm=norm, origin='lower')
plt.show()
In [5]:
mmirs = WFSFactory(wfs="mmirs")
In [6]:
mmirs_file = "/Users/tim/MMT/mmtwfs/notebooks/data.fits"
f9_file = "/Users/tim/MMT/wfsdat/20170318/TREY_p500_0000.fits"
ref_file = "/Users/tim/MMT/mmtwfs/mmtwfs/data/ref_images/mmirs_camera1_ref.fits"
#mmirs_file = "/Users/tim/MMT/wfsdat/20170414/rawdata/mmirs_wfs_0275.fits"
data = check_wfsdata(mmirs_file)
#cr_mask, data = detect_cosmics(data, sigclip=4., niter=6, cleantype='medmask')
In [8]:
implot(data)
In [252]:
d = data.copy()
sigma_clip = SigmaClip(sigma=3., iters=10)
bkg_estimator = MedianBackground()
bkg = Background2D(d, (12, 12), filter_size=(5, 5), bkg_estimator=bkg_estimator)
implot(d - bkg.background)
In [258]:
d = data.copy()
im = d - bkg.background
im[:5, :] = 0.0
im[:, :12] = 0.0
s = wfsfind(data, plot=True, fwhm=4.0, threshold=30.)
plt.show()
In [248]:
sf = s[s['flux'] > s['flux'].max()/3.]
sf.show_in_notebook()
Out[248]:
In [178]:
xspacing, yspacing = grid_spacing(im)
In [15]:
masked_data = data.copy()
masked_data[masked_data < 300.] = 0.0
xcen, ycen, plot = center_pupil(masked_data, mmirs.telescope.pupil_mask(), plot=True)
plt.show()
In [16]:
xcen, ycen
Out[16]:
In [259]:
apers = photutils.CircularAperture((s['xcentroid'], s['ycentroid']), r=(xspacing+yspacing)/4.)
In [260]:
masks = apers.to_mask(method='subpixel')
In [261]:
spot = np.zeros(masks[0].shape)
In [262]:
for m in masks:
subim = m.cutout(data)
try:
spot += subim
except:
pass
In [263]:
implot(spot)
In [264]:
g = Gaussian2D(amplitude=spot.max(), x_mean=spot.shape[1]/2, y_mean=spot.shape[0]/2) + Polynomial2D(degree=0)
fit_g = LevMarLSQFitter()
y, x = np.mgrid[:spot.shape[0], :spot.shape[1]]
In [265]:
pars = fit_g(g, x, y, spot)
In [266]:
pars.y_stddev_0.value
Out[266]:
In [189]:
pars
Out[189]:
In [190]:
ref = mmirs.modes['mmirs2']['reference']
In [191]:
refx = (xspacing/ref['xspacing'])*ref['apertures']['xcentroid']+305 #xcen
refy = (yspacing/ref['yspacing'])*ref['apertures']['ycentroid']+247 #ycen
ref_apers = photutils.CircularAperture(
(refx, refy),
r=(xspacing+yspacing)/4.,
)
In [192]:
ref_apers.plot()
#apers.plot()
implot(im)
plt.show()
In [193]:
src_coord = SkyCoord(x=s['xcentroid'], y=s['ycentroid'], z=0.0, representation='cartesian')
ref_coord = SkyCoord(x=refx, y=refy, z=0.0, representation='cartesian')
In [194]:
idx, sep, dist = match_coordinates_3d(src_coord, ref_coord)
In [195]:
src_coord.separation(ref_coord[idx]).std().to(u.deg)
Out[195]:
In [196]:
np.mean(abs(s['xcentroid'] - refx[idx]))
Out[196]:
In [197]:
np.mean(s['xcentroid'] - refx[idx])
Out[197]:
In [198]:
def fit_wfs_reference(pars, ref, spots):
xc = pars[0]
yc = pars[1]
model = SkyCoord(
representation='cartesian',
x=ref['xcentroid'] + xc,
y=ref['ycentroid'] + yc,
z=0.0
)
measured = SkyCoord(
representation='cartesian',
x=spots['xcentroid'],
y=spots['ycentroid'],
z=0.0
)
idx, ang, dist = match_coordinates_3d(measured, model)
stat = measured.separation(model[idx]).std().value
return stat
In [199]:
scale_ref = {
'xcentroid': (xspacing/ref['xspacing'])*ref['apertures']['xcentroid'],
'ycentroid': (yspacing/ref['yspacing'])*ref['apertures']['ycentroid']
}
pars = [xcen, ycen]
xopt = optimize.fmin_powell(fit_wfs_reference, pars, xtol=1e-4, args=(scale_ref, s))
xopt, xcen, ycen
Out[199]:
In [200]:
refx_match = refx[idx]
refy_match = refy[idx]
xoff = np.mean(s['xcentroid'] - refx_match)
yoff = np.mean(s['ycentroid'] - refy_match)
ref_coord = SkyCoord(x=refx+xoff, y=refy+yoff, z=0.0, representation='cartesian')
idx, sep, dist = match_coordinates_3d(src_coord, ref_coord)
xoff = np.mean(s['xcentroid'] - refx[idx])
yoff = np.mean(s['ycentroid'] - refy[idx])
In [201]:
dist
Out[201]:
In [202]:
ref_apers = photutils.CircularAperture(
(refx[idx]+xoff, refy[idx]+yoff),
r=(xspacing+yspacing)/4.,
)
plt.imshow(im, origin='lower')
ref_apers.plot()
#apers.plot()
plt.show()
In [203]:
len(ref['apertures'])
Out[203]:
In [204]:
ref_im = fits.open("/Users/tim/MMT/mmtwfs/mmtwfs/data/ref_images/mmirs_camera2_ref.fits")[0].data
In [221]:
ref_xsum = np.sum(ref_im, axis=0)
ref_xsum /= ref_xsum.max()
ref_ysum = np.sum(ref_im, axis=1)
ref_ysum /= ref_ysum.max()
xl = int(xcen - 2*xspacing)
xu = int(xcen + 2*xspacing)
yl = int(ycen - 2*yspacing)
yu = int(ycen + 2*yspacing)
xsum = np.sum(im[yl:yu, :], axis=0)
xsum /= xsum.max()
xtot = np.sum(im, axis=0)
xtot /= xtot.max()
ysum = np.sum(im[:, xl:xu], axis=1)
ysum /= ysum.max()
ytot = np.sum(im, axis=1)
ytot /= ytot.max()
xdiff = xtot - xsum
xdiff[xdiff < 0.25] = 0.0
ydiff = ytot - ysum
ydiff[ydiff < 0.25] = 0.0
In [222]:
#plt.plot(ysum)
plt.plot(xdiff)
plt.show()
In [207]:
ndimage.measurements.center_of_mass(xdiff)
Out[207]:
In [137]:
xcen, ycen
Out[137]:
In [223]:
detect_cosmics?
In [ ]: