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"]


WARNING:root:Unrecognized parameter: 'stpropflag'
WARNING:root:Unrecognized parameter: 'stpropflag'

In [194]:
pl.imshow(np.log(f[10], cmap="gray", interpolation="nearest")


Out[194]:
<matplotlib.image.AxesImage at 0x1247aac50>

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)


Removing existing file

In [199]:
check_call("bin/image2xy -O -w 1 -a 5 -p 8 tmp.fits", shell=True)


Out[199]:
0

In [200]:
coords = fitsio.read("tmp.xy.fits")

In [201]:
coords


Out[201]:
array([(12.758220672607422, 15.11616039276123, 19748.9453125, 503.1015625),
       (19.3499698638916, 3.028876781463623, 1229.046630859375, 503.5489501953125),
       (2.891984701156616, 19.11785316467285, 749.369384765625, 502.693115234375)], 
      dtype=[('X', '>f4'), ('Y', '>f4'), ('FLUX', '>f4'), ('BACKGROUND', '>f4')])

In [202]:
pl.imshow(img, cmap="gray", interpolation="nearest")
pl.plot(coords["X"]-1, coords["Y"]-1, "+r")


Out[202]:
[<matplotlib.lines.Line2D at 0x13fb1ea10>]

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


[[ 14.17405912   2.0545043   18.17858821]
 [ 14.15930798   2.04660051  18.16711595]
 [ 14.15753422   2.04567595  18.16566107]
 ..., 
 [ 14.06235694   2.02380664  18.06542857]
 [ 14.06132488   2.0172498   18.06373157]
 [ 14.06114889   2.01585051  18.06335276]]

Find centroid position for every star at every time


In [205]:
print np.shape(xvals)


(3753, 3)

In [208]:
n = 1000
pl.scatter(xvals[-n:, 0], yvals[-n:, 0], c=np.arange(n), lw=0, alpha=0.7)


Out[208]:
<matplotlib.collections.PathCollection at 0x13f9d5890>

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])


0.0
0.176963638257
0.0903595527393
0.176963638257
0.0
0.193970319172
0.0903595527393
0.193970319172
0.0

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]:
[<matplotlib.lines.Line2D at 0x13fd2cb90>]

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)


The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-282-da1b1ed30b72> in <module>()
      1 get_ipython().magic(u'load_ext autoreload')
      2 get_ipython().magic(u'autoreload 2')
----> 3 from kpsf._kpsf import (compute_model, run_photometry_all,
      4                         N_INT_TIME, N_PSF_COMP, N_STARS)

ImportError: cannot import name 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 [ ]: