In [1]:
%matplotlib inline
%config InlineBackend.figure_format = "retina"
from matplotlib import rcParams
rcParams["savefig.dpi"] = 100
rcParams["font.size"] = 20
In [169]:
import kplr
import fitsio
import numpy as np
import matplotlib.pyplot as pl
from subprocess import check_call
Load the pointing model
In [69]:
times = fitsio.read("data/pointing_model.fits", 0)
pointing = fitsio.read("data/pointing_model.fits", 1)
Load the data
In [193]:
# data = fitsio.read("data/ktwo202059521-c00_lpd-targ.fits.gz")
client = kplr.API()
tpf = client.k2_star(202135911).get_target_pixel_files()
data = tpf[0].read()
data_times = data["TIME"]
f = data["FLUX"]
ferr = data["FLUX_ERR"]
In [194]:
pl.imshow(np.log(f[10], cmap="gray", interpolation="nearest")
Out[194]:
Use image2xy to find the sources in this image
In [195]:
img = f[-10]
m = np.isfinite(img) * (img > 0.0)
img[~m] = np.median(img[m])
fitsio.write("tmp.fits", np.array(img, dtype=np.float64), clobber=True)
In [199]:
check_call("bin/image2xy -O -w 1 -a 5 -p 8 tmp.fits", shell=True)
Out[199]:
In [200]:
coords = fitsio.read("tmp.xy.fits")
In [201]:
coords
Out[201]:
In [202]:
pl.imshow(img, cmap="gray", interpolation="nearest")
pl.plot(coords["X"]-1, coords["Y"]-1, "+r")
Out[202]:
In [203]:
from kpsf.c3k import fit_5x5, fit_3x3
In [204]:
lt = np.shape(f)[0]
yvals = np.nan + np.zeros((lt, len(coords)))
xvals = yvals + 0.0
for t in xrange(lt):
img = f[t]
m = np.isfinite(img)
m[m] *= img[m] > 0.0
if not np.any(m):
continue
img[~m] = np.median(img[m])
for k, c in enumerate(coords):
yi = int(min(max(3, np.round(c["X"] - 1)), img.shape[1]-4))
xi = int(min(max(3, np.round(c["Y"] - 1)), img.shape[0]-4))
img_5x5 = img[xi-2:xi+3, yi-2:yi+3]
xi2, yi2 = np.unravel_index(np.argmax(img_5x5), img_5x5.shape)
xi2, yi2 = xi2 - 2, yi2 - 2
ox, oy = fit_3x3(img[xi + xi2 - 1:xi + xi2 + 2, yi+yi2-1:yi+yi2+2])
# pl.imshow(img.T, cmap="gray", interpolation="nearest")
# pl.plot(xi, yi, "+r")
# pl.plot(ox + xi + xi2, oy + yi + yi2, "+b")
xvals[t, k] = ox + xi + xi2
yvals[t, k] = oy + yi + yi2
print xvals
Find centroid position for every star at every time
In [205]:
print np.shape(xvals)
In [208]:
n = 1000
pl.scatter(xvals[-n:, 0], yvals[-n:, 0], c=np.arange(n), lw=0, alpha=0.7)
Out[208]:
In [251]:
from itertools import product
for i, j in product(range(3), range(3)):
d = yvals[-n:, i] - yvals[-n:, j]
print np.std(d[np.isfinite(d)])
pl.plot(data_times[-n:], yvals[-n:, i] - yvals[-n:, j])
In [210]:
dt = data_times[-n:]
inds = np.argmin(np.abs(times[:, None] - dt[None, :]), axis=0)
W = pointing[inds]
In [211]:
pl.plot(dt, W[:, 3])
Out[211]:
In [284]:
pred = np.nan + np.zeros((len(coords), n, 2))
a = None
mean_xs = np.zeros((len(coords), 2))
for c in range(len(coords)):
x = np.vstack((xvals[-n:, c], yvals[-n:, c])).T
m = np.all(np.isfinite(x), axis=1)
# The first time through, compute the A matrix.
if a is None:
a = np.linalg.solve(np.dot(W[m].T, W[m]), np.dot(W[m].T, x[m])).T
# Find the best fit x0 conditioned on this matrix.
a2 = np.array(a)
mean_xs[c] = np.mean(x[m] - np.dot(W[m, 1:], a[:, 1:].T), axis=0)
a2[:, 0] = mean_xs[c]
# Compute the prediction.
pred[c, m] = np.dot(W[m], a2.T)
pl.plot(data_times[-n:][m], pred[c, m, 1] - yvals[-n:, c][m])
In [278]:
pl.scatter(xvals[-n:, 1], yvals[-n:, 1], c=np.arange(n), lw=0, alpha=0.7)
#pl.scatter(pred[1, :, 0], pred[1, :, 1], c=np.arange(n), lw=0, alpha=0.7);
In [263]:
pl.figure(figsize=(8, 8))
for i in range(100):
pl.clf()
pl.imshow(np.log(f[-n + i].T), cmap="gray", interpolation="nearest")
pl.plot(pred[:, i, 0], pred[:, i, 1], "+r")
pl.savefig("frames/{0:05d}.png".format(i))
In [242]:
flux_images = np.ascontiguousarray(f[-n:], dtype=np.float64)
ferr_images = np.ascontiguousarray(ferr[-n:], dtype=np.float64)
In [282]:
%load_ext autoreload
%autoreload 2
from kpsf._kpsf import (compute_model, run_photometry_all,
N_INT_TIME, N_PSF_COMP, N_STARS)
In [294]:
import pickle
In [295]:
pickle.dump((mean_xs, a[:, 1:], W[:, 1:], times[-n:], f[-n:], ferr[-n:]), open("blugh.pkl", "wb"), -1)
In [ ]: