In [1]:
from __future__ import division, print_function
%matplotlib inline
%config InlineBackend.figure_format = "retina"
from matplotlib import rcParams
rcParams["savefig.dpi"] = 100
rcParams["font.size"] = 20
In [2]:
import h5py
import fitsio
import numpy as np
import matplotlib.pyplot as pl
from scipy.linalg import cho_solve, cho_factor
from astropy.io import fits
from astropy.wcs import WCS
First we'll read in the data:
In [3]:
tpf = fitsio.read("ktwo201466124-c01_lpd-targ.fits.gz")
Let's look at one of the images (the last one in the time series). It's in the FLUX column of the TPF file:
In [4]:
img = tpf["FLUX"][-1]
pl.imshow(img.T, cmap="gray", interpolation="nearest");
Here's some code for measuring the centroid using the 3x3 patch and modeling it as a 2D quadratic. I'll also calcuate the second derivatives of the flux in the 3x3 patch because this should be related to the shape of the PSF.
In [5]:
# These are some useful things to pre-compute and use later.
_x, _y = np.meshgrid(range(-1, 2), range(-1, 2), indexing="ij")
_x, _y = _x.flatten(), _y.flatten()
_AT = np.vstack((_x*_x, _y*_y, _x*_y, _x, _y, np.ones_like(_x)))
_ATA = np.dot(_AT, _AT.T)
factor = cho_factor(_ATA, overwrite_a=True)
# This function finds the centroid and second derivatives in a 3x3 patch.
def fit_3x3(img):
a, b, c, d, e, f = cho_solve(factor, np.dot(_AT, img.flatten()))
m = 1. / (4 * a * b - c*c)
x = (c * e - 2 * b * d) * m
y = (c * d - 2 * a * e) * m
dx2, dy2, dxdy = 2 * a, 2 * b, c
return [x, y, dx2, dy2, dxdy]
# This function finds the centroid in an image.
# You can provide an estimate of the centroid using WCS.
def find_centroid(img, init=None):
if init is None:
xi, yi = np.unravel_index(np.argmax(img), img.shape)
else:
xi, yi = map(int, map(np.round, init))
ox, oy = np.unravel_index(np.argmax(img[xi-1:xi+2, yi-1:yi+2]), (3, 3))
xi += ox - 1
yi += oy - 1
assert (xi >= 1 and xi < img.shape[0] - 1), "effed, x"
assert (yi >= 1 and yi < img.shape[1] - 1), "effed, y"
pos = fit_3x3(img[xi-1:xi+2, yi-1:yi+2])
pos[0] += xi
pos[1] += yi
return pos
We'll use astropy to get the WCS and estimate the location of the star.
In [6]:
with fits.open("ktwo201466124-c01_lpd-targ.fits.gz") as hdus:
hdr = hdus[2].header
wcs = WCS(hdr)
# The order is the opposite of what I normally use...
init = wcs.wcs_world2pix(hdr["RA_OBJ"], hdr["DEC_OBJ"], 0.0)[::-1]
Plot the initial centroid estimate:
In [7]:
img = tpf["FLUX"][-1]
pl.imshow(img.T, cmap="gray", interpolation="nearest")
pl.plot(init[0], init[1], "or");
And then refine it using the code above:
In [8]:
cx, cy = find_centroid(img, init)[:2]
In [9]:
img = tpf["FLUX"][-1]
pl.imshow(img.T, cmap="gray", interpolation="nearest")
pl.plot(init[0], init[1], "or")
pl.plot(cx, cy, ".b");
Now let's measure it at each time by looping over frames.
In [10]:
coords = np.nan + np.zeros((len(tpf), 5))
for i, img in enumerate(tpf["FLUX"]):
msk = np.isfinite(img)
if np.any(msk) and np.any(img[msk] > 0.0):
coords[i] = find_centroid(img, init=init)
In [11]:
print(coords.shape)
The elements of coords at each time are x, y, d^2f/dx^2, d^2f/dy^2, d^2f/dxdy.
Let's just use the photometry that I measured already...
In [12]:
data = fitsio.read("ktwo201466124-c01_lpd-lc.fits")
aps = fitsio.read("ktwo201466124-c01_lpd-lc.fits", 2)
Then we'll choose the "best" aperture size, mask the missing data and normalize the fluxes:
In [13]:
# Load the data.
time = data["time"]
flux = data["flux"][:, np.argmin(aps["cdpp6"])] # Lowest CDPP6.
qual = data["quality"]
# Mask missing/bad data.
good_mask = np.isfinite(time) & np.isfinite(flux) & (qual == 0)
good_time = time[good_mask]
good_flux = flux[good_mask] / np.median(flux[good_mask])
In [15]:
pl.plot(good_time, good_flux, "k")
pl.xlim(1980, 1990);
I think that what we want to do is build a non-linear model of coords (the centroids that we measured) that best predicts good_flux. Crossfield et al. use a Gaussian Process and I think that that would be a good idea but there are lots of other things to try.
Probably the best first bet is try plotting all the things against each other and see where we get...
Here's one to start with (color indicates time):
In [16]:
pl.scatter(coords[good_mask, 0], coords[good_mask, 1], c=good_time, edgecolor="none", alpha=0.5)
pl.xlabel("x [pix]")
pl.ylabel("y [pix]");
In [17]:
import george
from george.kernels import Matern32Kernel, ExpSquaredKernel
In [18]:
X = np.concatenate((coords[good_mask, :2], np.atleast_2d(good_time).T), axis=1)
X -= np.mean(X, axis=0)[None, :]
X /= np.std(X, axis=0)[None, :]
In [36]:
d = 1000.0 * np.ones(X.shape[1])
d[-1] = 100.0
kernel = np.var(good_flux) * ExpSquaredKernel(d, ndim=X.shape[1])
gp = george.GP(kernel, mean=1.0)
In [37]:
gp.compute(X, 1e-8)
In [38]:
mu = gp.predict(good_flux, X, mean_only=True)
In [39]:
mu.shape
Out[39]:
In [40]:
pl.plot(good_time, good_flux, "k")
pl.plot(good_time, mu, "r")
# pl.xlim(2040, 2050);
In [43]:
pl.plot(good_time, good_flux - mu, "k");
# pl.xlim(2040, 2050);
In [131]:
np.sum((X[500] - X[1000]) ** 2)
Out[131]:
In [ ]: