2d GP for LSST clouds

clouds.png

lsst_clouds.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import division, print_function

import numpy as np
import matplotlib.pyplot as pl

import george
from george import kernels

# Load the data.
data = np.load("cloudTest/cloudTest_03.npy")

# Select 1000 random points in the image as "stars".
K = 1000
xi = np.random.randint(data.shape[0], size=K)
yi = np.random.randint(data.shape[1], size=K)
X = np.vstack((xi, yi)).T
z = data[xi, yi]

# Set up the Gaussian Process.
kernel = np.var(z) * kernels.Matern32Kernel([1000., 1000.], ndim=2)
gp = george.GP(kernel, mean=np.mean(z))
gp.compute(X, 1e-8)

# Compute the prediction.
xx, yy = range(data.shape[0]), range(data.shape[1])
xgrid, ygrid = np.meshgrid(xx, yy, indexing="ij")
Xpred = np.vstack((xgrid.flatten(), ygrid.flatten())).T
Zpred = gp.predict(z, Xpred, mean_only=True).reshape(xgrid.shape)

# Plot the data.
fig = pl.figure(figsize=(10, 5))
ax = fig.add_subplot(131)
vmin, vmax = np.min(z), np.max(z)
ax.imshow(data, cmap="gray", interpolation="nearest", vmin=vmin, vmax=vmax)
ax.plot(xi, yi, ",r")
ax.set_xlim(0, data.shape[0])
ax.set_ylim(0, data.shape[1])
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.set_title("data")

# Plot the prediction.
ax = fig.add_subplot(132)
ax.imshow(Zpred, cmap="gray", interpolation="nearest", vmin=vmin, vmax=vmax)
ax.set_xlim(0, data.shape[0])
ax.set_ylim(0, data.shape[1])
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.set_title("prediction")

# Plot the residuals on the same scale.
ax = fig.add_subplot(133)
cax = ax.imshow(data - Zpred, cmap="gray", interpolation="nearest",
                vmin=vmin, vmax=vmax)
ax.set_xlim(0, data.shape[0])
ax.set_ylim(0, data.shape[1])
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.set_title("residuals")

pl.savefig("clouds.png")