Reprojection Demo#

A notebook on how to use the kbmod.reprojection module to reproject WorkUnits into a common WCS frame of reference.

[ ]:

from kbmod.work_unit import WorkUnit
from kbmod import reprojection

import matplotlib.pyplot as plt
from astropy.nddata import CCDData

import numpy as np
from os import path


# visualization utility.
def plot_images(obstimes, image1=None, image2=None, image3=None, image4=None):
    """Plot up to four images (`astropy.nddata.CCDData` objects) side by side"""

    min_val = -1.0
    fig = plt.figure(figsize=(10, 10))
    if image1:
        f1 = fig.add_subplot(2, 2, 1, projection=image1.wcs)
        f1.title.set_text(str(obstimes[0]))
        plt.imshow(image1.data, vmin=min_val, origin="lower")
        plt.grid(color="white", ls="solid")
    if image2:
        f2 = fig.add_subplot(2, 2, 2, projection=image2.wcs)
        f2.title.set_text(str(obstimes[1]))
        plt.imshow(image2.data, vmin=min_val, origin="lower")
        plt.grid(color="white", ls="solid")
    if image3:
        f3 = fig.add_subplot(2, 2, 3, projection=image3.wcs)
        f3.title.set_text(str(obstimes[2]))
        plt.imshow(image3.data, vmin=min_val, origin="lower")
        plt.grid(color="white", ls="solid")
    if image4:
        f4 = fig.add_subplot(2, 2, 4, projection=image4.wcs)
        f4.title.set_text(str(obstimes[3]))
        plt.imshow(image4.data, vmin=min_val, origin="lower")
        plt.grid(color="white", ls="solid")

Simple Reproject#

This the first and simplest type of reprojection that we do. Given a kbmod.search.WorkUnit with images of differing WCSs, choose a common WCS and reproject all the images to it.

First, we’ll need a WorkUnit. Let’s use the same one we usually use for reprojection testing.

[ ]:

project_root_dir = path.abspath((path.dirname(path.dirname(path.abspath("")))))
fake_data_loc = "tests/data/shifted_wcs_diff_dimms_tiled.fits"

wunit = WorkUnit.from_fits(path.join(project_root_dir, fake_data_loc))
wunit

Let’s take a peek into this WorkUnit and see what kind of data we have.

[ ]:

original_images = wunit.im_stack.get_images()

o_image0 = CCDData(original_images[0].get_science().image, unit="adu")
o_image0.wcs = wunit.get_wcs(0)

o_image1 = CCDData(original_images[1].get_science().image, unit="adu")
o_image1.wcs = wunit.get_wcs(1)

o_image2 = CCDData(original_images[2].get_science().image, unit="adu")
o_image2.wcs = wunit.get_wcs(2)

o_image3 = CCDData(original_images[3].get_science().image, unit="adu")
o_image3.wcs = wunit.get_wcs(3)

plot_images(wunit.get_all_obstimes(), o_image0, o_image1, o_image2, o_image3)

A couple of important attributes to point out:

  • Each images has a different WCS. The center ra/dec value shifted up and to the right ~5 pixels in each successive images, except for the last one which is below image 3.

  • The obstime is increasing for each one, except for the last one which has the same obstime as image 3.

  • They all have a synthetic object in them, moving across the field of view. The last image has a presumambly different object.

  • Image 3 and 4 are from the same “observation”, with the same obstime and aligned, to simulate the case where we have multiple tiled detectors from a given observation.

    • This is done mainly to handle the LSST case where observations will have to be built out of many different detectors as everything isn’t neatly aligned like in DEEP.

With that, let’s reproject the images and see how reprojection behaves.

Reprojection#

we’ll define a common_wcs, in this case the WCS of the first image in the stack. From there, we can run the original work unit through our reprojection module and get a new WorkUnit, where everything has been reprojected into the provided common_wcs.

[ ]:

common = wunit.get_wcs(0)

uwunit = reprojection.reproject_work_unit(wunit, common)

And now let’s take a look at the results.

[ ]:

images = uwunit.im_stack.get_images()

image0 = CCDData(images[0].get_science().image, unit="adu")
image0.wcs = uwunit.wcs

image1 = CCDData(images[1].get_science().image, unit="adu")
image1.wcs = uwunit.wcs

image2 = CCDData(images[2].get_science().image, unit="adu")
image2.wcs = uwunit.wcs

plot_images(uwunit.get_all_obstimes(), image0, image1, image2)

All of the images are now in the same WCS, which required some shifting - it’s hard to see in this view but along the left edge of the images except for the first have KB_NO_DATA for a lot of the pixels now. In addition, there are now only 3 images, as the last two were combined into one image as they had the same obstime. Also note that the first image has had no changes made to it, as we used its WCS, resulting in a no-op.

Important Note: We were only able to combine the last two images into one because they weren’t overlapping. Current constraints mean that if the two images have the same obstime but an overlapping footprint, an error will be thrown.

The reprojection will also change the variance and mask layers of the original LayeredImage. See below:

[ ]:

original_img = wunit.im_stack.get_single_image(2)
o_d = original_img.get_mask().image
original_image2_mask = CCDData(o_d, unit="adu")
original_image2_mask.wcs = wunit.get_wcs(2)

image2_mask = CCDData(images[2].get_mask().image, unit="adu")
image2_mask.wcs = uwunit.wcs

plot_images([uwunit.get_all_obstimes()[-1]] * 2, original_image2_mask, image2_mask)

Looking at the difference between the original mask and the new one, there was a line in the original mask that transfered over (the one around 02”). You can also see more clearly the parts with no data in the new image, on the left side as well as the line on the bottom that represents the “chip gap” between image 3 and 4 in the original ImageStack.

Barycentric Projection#

Coming soon!