Using the image cutout tool with DP0.2¶

Rubin Observatory logo, a graphical representation of turning stars into data.
Contact author: Christina Williams
Last verified to run: 2025-04-30
LSST Science Pipelines version: Weekly 2025_17
Container Size: medium
Targeted learning level: beginner

Description: This notebook demonstrates how to use the Rubin Image Cutout Service.

Skills: Run the Rubin Image Cutout Service for visual inspection of small cutouts of LSST images.

LSST Data Products: Images (deepCoadd, calexp), catalogs (objectTable, diaObject, truthTables, ivoa.ObsCore).

Packages: PyVO, lsst.rsp.get_tap_service, lsst.pipe.tasks.registerImage, lsst.afw.display

Credit: This notebook builds on technical exploration of the cutout service developed by Leanne Guy (in preparation), and builds on an earlier notebook written by Alex Drlica-Wagner linked here. This notebook additionally includes contributed functions by Jeff Carlin and Ryan Lau for time-variable sources.

Get Support: Find DP0-related documentation and resources at dp0.lsst.io. Questions are welcome as new topics in the Support - Data Preview 0 Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.

1. Introduction¶

This notebook will teach how to use the cutout service (which enables users to retrieve small image cutouts from LSST data) by demonstration of a few basic science applications. Since LSST images are quite large, in some cases it is desirable to be able to perform image cutouts on the server side to avoid transferring large amounts of data. This can speed up visual inspection and other analyses requiring images of individual objects or fields that may be much smaller than the patch and tract sizes that are created by the LSST camera and pipelines.

The International Virtual Observatory Alliance (IVOA) co-ordinates the community efforts of astronomical missions and archives to develop and maintain the Virtual Observatory (VO) standards. The VO standards enable interoperability between astronomical archives. IVOA also provides a Server-side Operations for Data Access (SODA), which is a protocol for remote data processing operations. This protocol allows users to perform computations (pixel operations, image transformations, etc) on the remote server, which avoids unnecssary data movement. The LSST architecture has a "VO-first" approach, meaning that VO standards are implemented in all applicable services, enabling the use of VO tools such as the image cutout service to access LSST data. The procedure is to identify the remote web location of the image of interest (called a datalink), and use a web service that creates a cutout from that linked data remotely, to transfer and save locally on the Rubin Science Platform.

Further details and information can be found at the IVOA data link documentation, where it says Access Data Services. Rubin-specific documentation for these can also be found in this document describing the RSP DataLink service implementation strategy.

1.1 Import packages¶

Import general python packages, LSST packages, PyVO packages, and Astropy packages.

In [1]:
import time
import numpy as np
import uuid
import os
import glob
import math
import pandas
import matplotlib.pyplot as plt
import tempfile
import getpass
from PIL import Image
from IPython.display import Image as dimg

import lsst.geom as geom
import lsst.resources
import lsst.afw.display as afwDisplay
from lsst.afw.image import Exposure, ExposureF
from lsst.pipe.tasks.registerImage import RegisterConfig, RegisterTask
from lsst.rsp import get_tap_service
from lsst.rsp.utils import get_access_token
from lsst.afw.fits import MemFileManager

import pyvo
from pyvo.dal.adhoc import DatalinkResults, SodaQuery

import astropy
from astropy import units as u
from astropy.coordinates import SkyCoord, Angle
from astropy.io import fits
from astropy.wcs import WCS

1.2 Define parameters and functions¶

1.2.1 Set plotting parameters¶

Set the backend for afwDisplay to matplotlib. Set parameters for plots.

In [2]:
afwDisplay.setDefaultBackend('matplotlib')

params = {'axes.labelsize': 18,
          'font.size': 18,
          'legend.fontsize': 12,
          'xtick.major.width': 2,
          'xtick.minor.width': 1,
          'xtick.major.size': 10,
          'xtick.minor.size': 4,
          'xtick.direction': 'in',
          'xtick.top': True,
          'lines.linewidth': 2,
          'axes.linewidth': 2,
          'axes.labelweight': 2,
          'axes.titleweight': 2,
          'ytick.major.width': 2,
          'ytick.minor.width': 1,
          'ytick.major.size': 10,
          'ytick.minor.size': 4,
          'ytick.direction': 'in',
          'ytick.right': True,
          'figure.figsize': [6, 6],
          'figure.facecolor': 'White'
          }
plt.style.use('tableau-colorblind10')

plt.rcParams.update(params)

When displaying pandas tables in the notebook, only display up to 20 rows.

In [3]:
pandas.set_option('display.max_rows', 20)

1.2.2 Set up the temporary directory¶

The Rubin cutout service allows the user to save cutouts as fits files locally on the Rubin Science Platform (RSP).

RSP users should save temporary files in the shared /scratch directory, in sub-directories named with their RSP username.

Get the username and ensure an appropriate /scratch sub-directory exists.

In [4]:
username = getpass.getuser()
print('username: ', username)

userdir = '/scratch/' + username
if not os.path.exists(userdir):
    os.makedirs(userdir)
    print('Created ', userdir)
else:
    print('Directory already existed: ', userdir)

username:  douglasleetucker
Directory already existed:  /scratch/douglasleetucker

Within that folder, create a sub-folder named dp02_13a_temp in which to save files created by this tutorial. At the end of the notebook, the last cell will clear the contents and remove this temporary folder.

In [5]:
tempdir = userdir + '/dp02_13a_temp'
if not os.path.exists(tempdir):
    os.makedirs(tempdir)
    print('Created ', tempdir)
else:
    print('Directory already existed: ', tempdir)

Created  /scratch/douglasleetucker/dp02_13a_temp

Delete the username and userdir, but keep tempdir to be used.

In [6]:
del username, userdir

1.2.3 Create functions¶

Create the function plot_image to enable easy plotting of files from the cutout tool.

In [7]:
def plot_image(exposure: ExposureF):
    """Plot and image using lsst.awf.image package

   Parameters
    ----------
    exposure : `ExposureF`
        the image to plot from file in LSST awf image exposure class format

   Returns
    -------
    image for notebook display
    """

    fig, ax = plt.subplots()
    display = afwDisplay.Display(frame=fig)
    display.scale('asinh', 'zscale')
    display.mtv(exposure.image)
    plt.show()

Create the function get_cutout to retrieve the cutout and hold it in memory as an ExposureF.

In [8]:
def read_cutout_mem(sq):
    """Read the cutout into memory

    Parameters
    ----------
    sq : 'dict'
        returned from SodaQuery.from_resource()

    Returns
    -------
    exposure : 'ExposureF'
        the cutout in exposureF format
    """

    cutout_bytes = sq.execute_stream().read()
    sq.raise_if_error()
    mem = MemFileManager(len(cutout_bytes))
    mem.setData(cutout_bytes, len(cutout_bytes))
    exposure = ExposureF(mem)

    return exposure

Create the function make_image_cutout.

This is a wrapper function that will call the cutout tool and create image cutouts:

  1. define the location on the sky
  2. query the TAP service for the specifications of the dataId
  3. retrieve the datalink URL associated with the data
  4. create a cutout instance from the query result and the Datalink Resource
  5. return the cutout as an ExposureF format image
In [9]:
def make_image_cutout(tap_service, ra, dec, dataId, cutout_size=0.01,
                      imtype=None, filename=None):
    """Wrapper function to generate a cutout using the cutout tool

   Parameters
    ----------
    tap_service : an instance of the TAP service
    ra, dec : 'float'
        the ra and dec of the cutout center
    dataId : 'dict'
        the dataId of the image to make a cutout from. The format
        must correspond to that provided for parameter 'imtype'
    cutout_size : 'float', optional
        edge length in degrees of the cutout
    imtype : 'string', optional
        string containing the type of LSST image to generate
        a cutout of (e.g. deepCoadd, calexp). If imtype=None,
        the function will assume a deepCoadd.
    filename : 'string', optional
        filename of the resulting cutout (which has fits format)

   Returns
    -------
    exposure : 'ExposureF'
        the cutout in exposureF format
    """

    spherePoint = geom.SpherePoint(ra*geom.degrees, dec*geom.degrees)

    if imtype == 'calexp':

        query = "SELECT access_format, access_url, dataproduct_subtype, " + \
            "lsst_visit, lsst_detector, lsst_band " + \
            "FROM ivoa.ObsCore WHERE dataproduct_type = 'image' " + \
            "AND obs_collection = 'LSST.DP02' " + \
            "AND dataproduct_subtype = 'lsst.calexp' " + \
            "AND lsst_visit = " + str(dataId["visit"]) + " " + \
            "AND lsst_detector = " + str(dataId["detector"])
        results = tap_service.search(query)

    else:
        # Find the tract and patch that contain this point
        tract = dataId["tract"]
        patch = dataId["patch"]

        # add optional default band if it is not contained in the dataId
        band = dataId["band"]

        query = "SELECT access_format, access_url, dataproduct_subtype, " + \
            "lsst_patch, lsst_tract, lsst_band " + \
            "FROM ivoa.ObsCore WHERE dataproduct_type = 'image' " + \
            "AND obs_collection = 'LSST.DP02' " + \
            "AND dataproduct_subtype = 'lsst.deepCoadd_calexp' " + \
            "AND lsst_tract = " + str(tract) + " " + \
            "AND lsst_patch = " + str(patch) + " " + \
            "AND lsst_band = " + "'" + str(band) + "' "
        results = tap_service.search(query)

    # Get datalink
    dataLinkUrl = results[0].getdataurl()
    auth_session = service._session
    dl = DatalinkResults.from_result_url(dataLinkUrl,
                                         session=auth_session)

    # from_resource: creates a instance from
    # a number of records and a Datalink Resource.
    sq = SodaQuery.from_resource(dl,
                                 dl.get_adhocservice_by_id("cutout-sync"),
                                 session=auth_session)

    sq.circle = (spherePoint.getRa().asDegrees() * u.deg,
                 spherePoint.getDec().asDegrees() * u.deg,
                 cutout_size * u.deg)

    exposure = read_cutout_mem(sq)

    # cutout_bytes = sq.execute_stream().read()
    # mem = MemFileManager(len(cutout_bytes))
    # mem.setData(cutout_bytes, len(cutout_bytes))
    # exposure = ExposureF(mem)
    
    return exposure

Create three functions to warp and scale a calexp cutout (which can have any rotated orientation on the sky) and combine them into an animated gif, to be used in Section 4.

  • warp_image: rotate and scale the image to a common grid
  • get_minmax_xy: determine the bounding box of the cutouts
  • make_gif: combine warped images into an animated gif
In [10]:
def warp_img(ref_img, img_to_warp, ref_wcs, wcs_to_warp):
    """Warp and rotate an image onto the coordinate system of another image

   Parameters
    ----------
    ref_img: 'ExposureF'
        is the reference image for the re-projection
    img_to_warp: 'ExposureF'
        the image to rotate and warp onto the reference image's wcs
    ref_wcs: 'wcs object'
        the wcs of the reference image (i.e. ref_img.getWcs() )
    wcs_to_warp: 'wcs object'
        the wcs of the image to warp (i.e. img_to_warp.getWcs() )
   Returns
    -------
    warpedExp: 'ExposureF'
        a reprojected, rotated image that is aligned and matched to ref_image
     """
    config = RegisterConfig()
    task = RegisterTask(name="register", config=config)
    warpedExp = task.warpExposure(img_to_warp, wcs_to_warp, ref_wcs,
                                  ref_img.getBBox())

    return warpedExp
In [11]:
def get_minmax_xy(img, cutout_size):
    """Get the pixel dimensions of an image

   Parameters
    ----------
    img: 'ExposureF'
        is the input image to return the pixel coordinates of
    cutout_size: 'int'
        is the edge size of the image in pixels
   Returns
    -------
    minx, maxx, miny, maxy: 'int'
        the min and max x and y pixel coordinates for the input image
     """

    cutout_size = int(cutout_size)

    height = img.height
    width = img.width

    ceny = (height - 1) / 2
    cenx = (width - 1) / 2

    minx = int(cenx - ((cutout_size - 1) / 2))
    maxx = int(cenx + ((cutout_size - 1) / 2))
    miny = int(ceny - ((cutout_size - 1) / 2))
    maxy = int(ceny + ((cutout_size - 1) / 2))

    return {'minx': minx, 'maxx': maxx, 'miny': miny, 'maxy': maxy}
In [12]:
def make_gif(folder):
    """Generate a GIF for all *.png images contained in a folder

   Parameters
    ----------
    # folder: 'string'
        string containing the path to the folder
        default filename is animation.gif

   Returns
    -------
     """
    frames = [Image.open(img) for img in sorted(glob.glob(f"{folder}/*.png"))]
    frame_one = frames[0]
    frame_one.save(folder+"/animation.gif", format="GIF",
                   append_images=frames, save_all=True, duration=500, loop=1)

2. A step by step demonstration of how to use the Rubin Image Cutout Service¶

This section will demonstrate a simple example: how to create a cutout for a single part of sky from a deepCoadd.

2.1 Initiate the TAP service, and define sky coordinates for the image cutout¶

The TAP service is used to query the ivoa.Obscore table for the datalink (a web URL identifying where the data is hosted).

In [13]:
service = get_tap_service("tap")

First, define a point on the sky as the center of the image cutout. This example uses the galaxy cluster from DP0.2 Notebook Tutorial 03a. Once the RA and Dec are defined, create a SpherePoint object to define the location on the sky. Use these in a TAP query to identify the overlapping deepCoadd image.

In [14]:
ra = 55.7467
dec = -32.2862
coord = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, frame='icrs')
radius = .1 * u.deg

spherePoint = lsst.geom.SpherePoint(ra*geom.degrees, dec*geom.degrees)

The next cell shows the TAP query for the metadata that is associated with the image's remote location in the LSST data archive. The DP0.2 has a schema (table collection) called "ivoa", which contains a table called ivoa.ObsCore. The IVOA-defined obscore table contains generic metadata for the DP0.2 datasets. The table is accessible via ADQL queries via a TAP endpoint. The mechanism for locating the images LSST is to use the TAP service to query the ObsCore schema.

In [15]:
query = "SELECT lsst_patch, lsst_tract, s_region, " + \
        "access_format, access_url, dataproduct_subtype " + \
        "FROM ivoa.ObsCore " + \
        "WHERE dataproduct_type = 'image' " + \
        "AND obs_collection = 'LSST.DP02' " + \
        "AND dataproduct_subtype = 'lsst.deepCoadd_calexp' " + \
        "AND lsst_band = 'i' " + \
        "AND CONTAINS(POINT('ICRS', " + str(coord.ra.value) + \
        ", " + str(coord.dec.value) + "), " + \
        "s_region) = 1"
In [16]:
job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)
assert job.phase == 'COMPLETED'

Job phase is COMPLETED
In [17]:
results = job.fetch_result()
In [18]:
results.to_table()
Out[18]:
Table length=1
lsst_patchlsst_tracts_regionaccess_formataccess_urldataproduct_subtype
int64int64objectobjectobjectobject
174431POLYGON ICRS 55.514085 -32.322253 55.790197 -32.322253 55.789845 -32.088924 55.514437 -32.088924application/x-votable+xml;content=datalinkhttps://data.lsst.cloud/api/datalink/links?ID=butler%3A//dp02/20d28216-534a-4102-b8a7-1c7f32a9b78clsst.deepCoadd_calexp

In the above table, the access_url contains the web URL datalink for the requested deepCoadd. This datalink will be needed to generate the image cutout.

The identified Patch and Tract will be used to define the unique dataId for that location, and once a filter (band) is included, this defines a unique deepCoadd in the LSST image database. Below, construct the dataId.

In [19]:
tract = results['lsst_tract'][0]
patch = results['lsst_patch'][0]

dataId = {'band': 'i', 'tract': tract,
          'patch': patch}

2.2 Generating an image cutout¶

First, extract the datalink (dataLinkUrl) from the query result using the getdataurl method. Next, extract the session authentication, which is needed for reuse of the datalink. Finally, create a DatalinkResults object to be able to access this URL, which will be stored as dl_results and available for approximately 15 minutes, in a format that can be used by the IVOA tools below.

In [20]:
dataLinkUrl = results[0].getdataurl()

auth_session = service._session
dl_results = DatalinkResults.from_result_url(dataLinkUrl,
                                             session=auth_session)

f"Datalink status: {dl_results.status}. Datalink service url: {dataLinkUrl}"
Out[20]:
"Datalink status: ('OK', 'QUERY_STATUS not specified'). Datalink service url: https://data.lsst.cloud/api/datalink/links?ID=butler%3A//dp02/20d28216-534a-4102-b8a7-1c7f32a9b78c"

Lastly, call the Rubin Image Cutout Service. It is done by feeding the data link created above (called dl_results) to following function:

In [21]:
sq = SodaQuery.from_resource(dl_results,
                             dl_results.get_adhocservice_by_id("cutout-sync"),
                             session=auth_session)

The variable sq now holds the result of the SODA query using the data link (which currently still points the full LSST patch and tract deepCoadd image, at its remote location in the database). The cell below will now demonstrate how to extract a cutout from sq.

2.2.1 Defining a square cutout using a subtended circle¶

Only two shape definitions are supported: a circle function, and a polygon function can be used to define the cutout dimensions. These shape definitions do not produce circle or polygon cutouts, but rather are methods for defining the edges of cutouts with 4 sides. In the case of circle, the resulting cutout is always a square, with edge size that is the same as the circle diameter. In the case of a polygon, either a square or a rectangular cutout will result, depending on whether the length and width edge dimensions are different values. Only cutouts with 4 corners and 90 degree angles are supported.

In [22]:
sphereRadius = 0.05 * u.deg
sq.circle = (spherePoint.getRa().asDegrees() * u.deg,
             spherePoint.getDec().asDegrees() * u.deg,
             sphereRadius)

2.2.2 Retrieve and plot the cutout¶

The cutout procedure is performed remotely.

Bytes are returned, and the user has two options:

  1. hold the bytes in memory and plot the image, or
  2. save the bytes to a file, and then read the file and plot the image.

Option 1: retrieve the bytes to memory, convert to ExposureF format, and display the cutout. This is what the function read_cutout_mem does.

In [23]:
help(read_cutout_mem)

Help on function read_cutout_mem in module __main__:

read_cutout_mem(sq)
    Read the cutout into memory

    Parameters
    ----------
    sq : 'dict'
        returned from SodaQuery.from_resource()

    Returns
    -------
    exposure : 'ExposureF'
        the cutout in exposureF format
In [24]:
exposure = read_cutout_mem(sq)
plot_image(exposure)
No description has been provided for this image

Figure 1: The deepCoadd cutout image, displayed in grayscale with a scale bar at right, and image axes in pixels. The field is centered on a rich galaxy cluster.

Clear the memory.

In [25]:
del exposure

Option 2: save the cutout to a file named "cutout-circle.fits" in the temporary directory. Then load it as an ExposureF and plot it.

In [26]:
sodaCutout = os.path.join(tempdir, 'cutout-circle.fits')
response_stream = sq.execute_stream()
sq.raise_if_error()
with open(sodaCutout, 'bw') as f:
    f.write(response_stream.read())
plot_image(ExposureF(sodaCutout))
No description has been provided for this image

Above is identical to Figure 1.

Future planned options for the Rubin cutout service, including the potential to retrieve other image formats such as jpeg, are listed at the Rubin Science Platform image cutout implementation strategy document.

Option "2b": If the file will not be needed again, save it as a temporary file which is automatically deleted. Uncomment and execute the following cell and see the exact same image is displayed as in Figure 1.

In [27]:
# with tempfile.NamedTemporaryFile(dir=tempdir, suffix=".fits") as temp:
#     response_stream = sq.execute_stream()
#     sq.raise_if_error()
#     temp.write(response_stream.read())
#     temp.flush()
#     plot_image(ExposureF(temp.name))

2.2.3 Using polygon to define the image cutout shape instead of a circle¶

It is also possible to define the cutout geometry using a polygon, which enables the cutout to be rectangular, but not necessarily be square. For this, use polygon, which takes as input the four corners in celestial coordinates. A minimum of 3 vertices are required (the line from the last vertex back to the first is implicit) Vertices must be ordered in the counter-clockwise direction. For example: a polygon is defined as a set of 4 (x,y) coordinates from (12,34) to (14,34) to (14,36) to (12,36) and (implicitly) back to (12,34) as:

POLYGON=12 34 14 34 14 36 12 36

Since the center of the galaxy cluster is already defined in RA and Dec in the cells above (spherePoint), this example will define each x,y set as RA+/-sphereRadius and Dec+/-sphereRadius.

In [28]:
sq2 = SodaQuery.from_resource(dl_results,
                              dl_results.get_adhocservice_by_id("cutout-sync"),
                              session=auth_session)
sphereRadius1 = 0.03 * u.deg
sphereRadius2 = 0.01 * u.deg

sq2.polygon = (spherePoint.getRa().asDegrees() * u.deg - sphereRadius1,
               spherePoint.getDec().asDegrees() * u.deg - sphereRadius2,
               spherePoint.getRa().asDegrees() * u.deg - sphereRadius1,
               spherePoint.getDec().asDegrees() * u.deg + sphereRadius2,
               spherePoint.getRa().asDegrees() * u.deg + sphereRadius1,
               spherePoint.getDec().asDegrees() * u.deg + sphereRadius2)

exposure = read_cutout_mem(sq2)
plot_image(exposure)
No description has been provided for this image

Figure 2: A rectangular polygon cutout of the image displayed in Figure 1.

In [29]:
del exposure

There is an important difference to note between the circle and polygon shape definitions. The angular distance on the sky that defines the circular cutout size already accounts for the difference in angular distance in the RA direction is smaller by a factor of cos(declination), where declination is in units radians. The difference increases with higher declination. However, the polygon definition does not automatically account for this cosine factor. Thus, circle and polygon cutout definitions using the same cutout edge length will not match size in the RA direction (for deepCoadds, the x-axis). The cell below demonstrates how to make this correction to the polygon cutout definition to create symmetric cutouts with polygon. Here sphereRadius comes from the circle definition above.

In [30]:
a = Angle(spherePoint.getDec().asDegrees(), u.deg)
cosd = np.cos(a.radian)

sq2.polygon = (spherePoint.getRa().asDegrees() * u.deg - sphereRadius/cosd,
               spherePoint.getDec().asDegrees() * u.deg - sphereRadius,
               spherePoint.getRa().asDegrees() * u.deg - sphereRadius/cosd,
               spherePoint.getDec().asDegrees() * u.deg + sphereRadius,
               spherePoint.getRa().asDegrees() * u.deg + sphereRadius/cosd,
               spherePoint.getDec().asDegrees() * u.deg + sphereRadius,
               spherePoint.getRa().asDegrees() * u.deg + sphereRadius/cosd,
               spherePoint.getDec().asDegrees() * u.deg - sphereRadius)

exposure = read_cutout_mem(sq2)
plot_image(exposure)
No description has been provided for this image

Figure 3: The same cutout as displayed in Figure 1, but using the polygon functionality and accounting for the $cos({\rm dec})$ term.

In [31]:
del exposure

2.3 Test out the cutout wrapper function "make_image_cutout"¶

All of the above steps from section 2.2 have been compiled in a wrapper function called make_image_cutout that is defined in Section 1. In the rest of this notebook, cutouts will be generated using this function for simplicity. As defined, it assumes the circular cutout definition demonstrated above. Thus, the function requires as input the TAP service, the center ra/dec of the cutout, the dataId and imtype (Section 4 will demonstrate how to do this for calexp, not just deepCoadd) and the size of the cutout (i.e. what was defined as sphereRadius above). The next cell demonstrates how to run all the steps by calling make_image_cutout and plotting the result.

In [32]:
imtype = 'deepCoadd'
test = make_image_cutout(service, ra, dec, dataId=dataId,
                         imtype=imtype, cutout_size=0.005)
plot_image(test)
No description has been provided for this image

Figure 4: A small image cutout, created using the make_image_cutout function.

3. Science application: validating a galaxy sample based on color-color selection¶

This section will demonstrate a simple example of a science use-case for the Rubin Image Cutout Service.

Lyman-break color selections are a common way to identify galaxies at high redshifts. The selections make use of the absorption by intergalactic hydrogen in the foreground of all galaxy light blueward of the Lyman limit (912A) or at very high redshifts, Lyman-alpha (1216A). More details about why this works can be found at this blog-post. DP0.2 includes mock galaxies up to redshift (z) of 3, where the Lyman-break falls between the LSST u and g bands. Such galaxies can be identified by their very red u-g colors, and since intergalactic hydrogen absorbs all the light blueward of the Lyman break, true z~3 galaxies should be undetected in the LSST u-band. The cell below will perform a search for redshift of 3 galaxies as u-band dropouts, and will create cutouts in each of the LSST filters in order to visually inspect the Lyman-break galaxy (LBG) candidates to confirm the real candidates show no u-band flux.

3.1 The u-band dropout LBG color selection¶

To search for $z\sim3$ LBGs in DP0.2, first define a u-band color selection based on the existing literature that uses filter sets similar to that of LSST. Below outlines the set of colors used in Steidel et al. 2003 (see also Boutsia et al 2014, among other references. Many selections exist, where changes to colors would alter the redshift selection function, or may differ due to differences in the broad band filter shapes of other instruments).

U − G >= (G − R) + 1.0

G − R <= 1.2

where capital letters U, G, and R indicate the total AB magnitudes in the LSST filters u, g, r. Since DP0.2 is only 5-year depth mock data in the LSST filters, our search will also include a g-band magnitude that is relatively bright (G < 24.5 ABmag) in order to obtain high-quality candidates.

High redshift galaxies are typically small in size. For this exercise, use aperture photometry with a relatively small aperture diameter (9 pixels) in order to obtain high S/N measurements and exclude neighbors.

3.1.1 Amending a typical search because of the rarity of z~3 galaxies in the mock data.¶

The mock galaxies that go into DP0.2 are created using the LSST Catalog Simulator (catSim), and are based on an empirical model of the evolution of galaxies outlined in Korytov et al. 2019. The input model is based on realistic number densities, fluxes, and redshift distributions across cosmic time, and utilizes the UniverseMachine methodology to assign empirical properties (see Behroozi et al. 2019).

Typically at any redshift, bright galaxies are rarer than faint galaxies, which are much more common. And typically in any patch of the sky, nearby galaxies will be much more common above the detection limit of a survey than distant $z\sim3$ galaxies. This means that $z\sim3$ galaxies that are bright enough to meet the criteria we want (G$ < 24.5$) will require a large search area to identify a statistical sample. So, in this example, we will first identify a parent sample of bright (G < 24.5) galaxies over a very large area (4 degrees) to ensure the parent sample contains enough rare distant galaxies.

In [33]:
center_coords = "62, -37"
ralbg = 62.
declbg = -37.
radius = "4.0"
max_rec = 500000

query = "SELECT TOP " + str(max_rec) + " " + \
        "objectId, coord_ra, coord_dec, detect_isPrimary, patch, tract, " + \
        "u_ap09Flux, u_ap09FluxErr, " + \
        "scisql_nanojanskyToAbMag(u_ap09Flux) as umag, " + \
        "scisql_nanojanskyToAbMag(r_ap09Flux) as rmag, " + \
        "scisql_nanojanskyToAbMag(g_ap09Flux) as gmag, " + \
        "scisql_nanojanskyToAbMag(u_ap09FluxErr) as umagErr " + \
        "FROM dp02_dc2_catalogs.Object " + \
        "WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), " + \
        "CIRCLE('ICRS', " + center_coords + ", " + radius + ")) = 1 " + \
        "AND detect_isPrimary = 1 AND detect_fromBlend = 0 " + \
        "AND g_ap09Flux/g_ap09FluxErr > 10 AND " + \
        "r_ap09Flux/r_ap09FluxErr > 10" + \
        "AND scisql_nanojanskyToAbMag(g_ap09Flux) < 24.5 "
print(query)

SELECT TOP 500000 objectId, coord_ra, coord_dec, detect_isPrimary, patch, tract, u_ap09Flux, u_ap09FluxErr, scisql_nanojanskyToAbMag(u_ap09Flux) as umag, scisql_nanojanskyToAbMag(r_ap09Flux) as rmag, scisql_nanojanskyToAbMag(g_ap09Flux) as gmag, scisql_nanojanskyToAbMag(u_ap09FluxErr) as umagErr FROM dp02_dc2_catalogs.Object WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), CIRCLE('ICRS', 62, -37, 4.0)) = 1 AND detect_isPrimary = 1 AND detect_fromBlend = 0 AND g_ap09Flux/g_ap09FluxErr > 10 AND r_ap09Flux/r_ap09FluxErr > 10AND scisql_nanojanskyToAbMag(g_ap09Flux) < 24.5
In [34]:
job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)

Job phase is COMPLETED
In [35]:
results = job.fetch_result()
In [36]:
results.to_table()
Out[36]:
Table length=162452
objectIdcoord_racoord_decdetect_isPrimarypatchtractu_ap09Fluxu_ap09FluxErrumagrmaggmagumagErr
degdegnJynJy
int64float64float64boolint64int64float64float64float64float64float64float64
140584436336742478260.4766515-40.3873187True373261626.996277365.176334324.40683759430625323.19071463739498623.93025507983664526.86477517281199
140584436336742480360.5614135-40.3870738True373261537.226964464.524202624.5746054935626723.6467172961431224.25461928683714826.875693384572
140584436336742453460.6323949-40.3923925True373261531.822775962.801712224.58558266865873424.4949472374579524.35282832942051526.905071289189046
140584436336742459160.5135566-40.3911427True3732612743.732980469.188872722.80414539092406322.7615161457777622.8540547016501526.79990936333654
140584436336742458860.6161867-40.3912253True3732611329.303045668.077458223.59094000056114423.9716759436843123.84872004084805426.81749166955609
140584436336742457560.6735168-40.3911105True3732611078.470027364.909701623.8179798006249622.33814027653878823.2236054658406926.86922596845902
140584436336742516260.6904022-40.3791836True373261686.576417968.159274824.30827779382746524.47780617001177824.4310375830596326.816187597541248
140584436336742515960.5817868-40.3794797True373261123.17201763.354507526.1737198668507522.67464516244155324.0842684987559126.895556202040645
140584436336742418460.6197004-40.3992424True3732611877.466747164.354402523.21606936574139721.25647875834772222.2780621766218926.87855434388595
140584436336742408860.5755128-40.401769True373261168.673463664.639110825.83238309255788322.84369983527811624.11645047911073226.873761566277885
140584436336742424060.7146587-40.3980953True3732612145.354821467.613634323.07125217334582623.3359050299501423.1299578877615826.82491429919181
140584436336742607360.6020345-40.3598672True373261627.359461463.441077924.406208869734424.2496390517148324.35200366718582326.89407361679371
140584436336742616160.5116247-40.3584183True373261275.348906764.572701325.3002916089709123.87002492009955324.43018513972059226.87487761292512
140584436336742614560.7217293-40.3578056True3732611244.666219167.851147223.6623677431838222.12631521679002623.14040550893273726.821107012662864
140584436336742604660.7187666-40.3600953True373261963.393344770.157274523.94049089549523424.0756181434690824.10061032361801426.784818227433217
140584436336742598060.7174846-40.3616604True3732613432.748461272.121921922.56089504715594522.7908727278852122.64303006603330526.754831771865092
140584436336742634660.5295535-40.3545276True373261539.420776466.366835224.57018082554715523.72400262110858224.2803712944167226.84512222958584
140584436336742538560.5289979-40.3749508True373261850.146960963.355426724.0762649833874924.24212369491348824.1198585705480126.895540449389117
140584436336742537160.6617351-40.3748337True373261635.626701967.399218224.3919946672541324.48188842389954424.45236645390643226.828362852630363
140584436336742545260.5290511-40.3733029True373261907.936589564.522803224.0048612041154624.19637501357378524.13833531318544326.875716932258914
....................................
182386988815483662663.4108633-33.2069286True3142312219.619113972.802439423.034303859958323.26159897505022323.12286279763072426.74463517121569
182386988815483646863.3998959-33.2101797True3142311133.176816669.959054323.76425579797537323.4281327572369123.70217305493951826.7878901747654
182386988815483640663.5156069-33.2111139True314231960.416749669.463428423.9438506862124.1185835908088924.0996149973357526.795609464222455
182386988815483651863.5514669-33.2092385True314231595.681003877.886749524.46246562320588724.24316903392124.43465742308871626.671341051213567
182386988815483543063.4631691-33.2300367True314231276.647490368.238546425.29518316319808223.58501759028672524.46414605086945626.81492558242895
182386988815483537863.5823381-33.2311943True3142311019.161808268.69524323.87939214835081720.61183811119086721.75917555028383526.80768333967935
182386988815483503063.4083057-33.2378641True314231924.911852673.860011523.9847491377566323.99487589965179523.92166922154874726.728976572663512
182386988815483502563.604175-33.2378321True314231654.161768968.685056224.36078710222790322.75986087109916823.82137578456538326.807844355161997
182386988815483575463.4656919-33.2242833True3142312132.313559170.599901723.07787232864708223.36587977245224223.1483184374525726.77798975909773
182386988815483554963.5154826-33.2275442True3142311215.117914870.311269623.68845394042712723.30920775146770323.5333922877026626.78243765000957
182386988815483557163.5450485-33.2274474True3142311004.50218766.605452623.89512278342422722.69534409289391523.01476188604932826.841225540544187
182386988815483555263.535231-33.2277942True314231531.669720567.356955124.58589518198783423.28766601294856324.1256685634982526.82904388383177
182386988815483564963.402924-33.2259625True3142311283.764103572.591926523.62878693046605623.6294872492198623.5208226299712226.74777919451258
182386988815483564463.4791333-33.226106True3142312646.058952573.741549922.8435012105940523.09312520249767823.06858412071090226.73071934657985
182386988815483561563.5617912-33.2266002True314231688.281319967.948590424.30558504282509324.2238431767697624.33195478882182526.819548870796396
182386988815483772663.4722067-33.1856639True314231830.995172471.401638324.10100374801905724.13924349572573324.0875587020655926.765729558217583
182386988815483756963.4097775-33.1882176True314231997.801884573.599161523.9023892004128723.55478190512408423.6022922846081826.73281783365345
182386988815483709663.4147072-33.1983661True314231236.041689369.047643525.46752821467246321.88745587434656423.08553175048035326.802127846747055
182385229596878574063.9396161-33.3378239True294231187.58735570.50556625.7169911003776223.77685184302609424.3134629533298426.77944149164673
182385229596878575163.9433615-33.3376814True294231526.937790470.697852524.59560163502722624.1803412081624824.46039317390031426.776484445053885
182385229596878506064.0015246-33.3508485True2942311144.890684369.02369623.75308994591224423.42932125032839323.52264958815365526.802504473324568

The next cells create a plot of the U-G vs G-R colors of the parent sample, and identify the high-redshift galaxy candidates based on the color selection discussed above. For sources with low S/N fluxes in u-band, we cannot robustly recover the U-G color, and instead we get a lower limit on how red it is (owing to the upper limit to the u-band flux). Thus, the U-G color can be overestimated if the u-band flux is not significant. Thus, sources with low S/N u-band flux will be reset to an error floor that corresponds the noise level of the u-band imaging so that it's color is not overestimated. To do this, define a new parameter called umag to store the u-band magnitude for assessing colors. In this case, if sources have signal to noise less than 1 in U_band, umag will instead hold the umagErr from the object table).

In [37]:
umag = np.zeros(len(results['umag']), dtype='float')
S2N = np.asarray(results['u_ap09Flux']/results['u_ap09FluxErr'], dtype='float')

for i in range(len(umag)):
    if S2N[i] < 1:
        umag[i] = results['umagErr'][i]
    else:
        umag[i] = results['umag'][i]

Further, this simple example uses a restricted color selection that is more narrow and blue in G-R than typical selections (G-R$<$0.2; this may exclude real z=3 galaxies that are redder in G-R, for example galaxies that have older stars or are dusty; but the restriction is helpful for this simple example because it lowers the number of interloper sources). The LBG candidates (indexed from the array as whlbg) identified using the color selection (black dashed lines) are in orange, with the full query of all objects in blue.

In [38]:
GRlim = 0.2
GRfloor = 1.2

whlbg = np.where(((umag-results['gmag'])
                  >= ((results['gmag'] - results['rmag']) + GRfloor))
                 & ((results['gmag'] - results['rmag']) < GRlim)
                 & ((results['gmag'] - results['rmag']) > 0))[0]

GmR = np.arange(0, GRlim + 0.01, 0.1)
plt.plot(GmR, GmR + GRfloor, '--', color='k')
plt.plot([0, 0], [GRfloor, 3], '--', color='k')
plt.plot([GRlim, GRlim], [GRlim + GRfloor, 3], '--', color='k')

GminusR = results['gmag']-results['rmag']
UminusG = umag-results['gmag']

plt.plot(GminusR, UminusG, '.', label='all objects', alpha=.1)
plt.plot(GminusR[whlbg], UminusG[whlbg], '.',
         label='color-selected LBG candidates')

plt.xlabel('G - R')
plt.ylabel('U - G')
plt.legend()
plt.show()
No description has been provided for this image

Figure 5: The $g-r$ versus $u-g$ color of all objects (blue) and color-selected LBG candidates (orange), with the box defining LBG candidates marked with black dashed lines.

Some fraction of the orange LBG candidates might be interlopers, entering the selection box for various reasons (e.g. photometric scatter). The next two cells below will plot image cutouts of each LBG candidate for visual inspection. Since this is a simulated dataset, the first cell performs query3 to the TAP service to pull the intrinsic redshift ts_redshift for each object from the truth table. This will be used to compare with the expected redshift $z\sim3$ of selected objects in the image cutout captions (see DP0.2 Tutorial Notebook 08 for more information on truth tables).

To limit the compute time (dominated by creating the cutouts), select just the first 5 galaxy candidates among the sample. This number can be altered using n_plot.

In [39]:
n_plot = 5

obj_id = str(tuple(int(objid) for objid in results['objectId'][whlbg][0:n_plot]))

query3 = "SELECT mt.id_truth_type AS mt_id_truth_type, " + \
         "obj.objectId AS obj_objectId,  " + \
         "ts.redshift AS ts_redshift " + \
         "FROM dp02_dc2_catalogs.MatchesTruth AS mt " + \
         "JOIN dp02_dc2_catalogs.TruthSummary AS ts " + \
         "ON mt.id_truth_type=ts.id_truth_type " + \
         "JOIN dp02_dc2_catalogs.Object AS obj " + \
         "ON mt.match_objectId=obj.objectId " + \
         "WHERE obj.objectId IN "+obj_id+" " + \
         "AND ts.truth_type=1 " + \
         "AND obj.detect_isPrimary=1 "

results3 = service.search(query3)
results3.to_table()
Out[39]:
Table length=5
mt_id_truth_typeobj_objectIdts_redshift
objectint64float32
9762571970_114865661090322262621.98357
9706148291_114862846340555160572.83135
8769499238_114851411419626224182.9079
8769504807_115668040695808055002.91797
7957723427_115675957179528017042.97819

Perform an inner join between results3 with results using the obj_objectId and objectId columns to ensure the correct rows from results3 are matched with the correct rows from results in the rest of this section.

In [40]:
joined_results = astropy.table.join(results3.to_table(), results.to_table(), 
                                    keys_left='obj_objectId',
                                    keys_right='objectId',
                                    join_type='inner')

joined_results
Out[40]:
Table length=5
mt_id_truth_typeobj_objectIdts_redshiftobjectIdcoord_racoord_decdetect_isPrimarypatchtractu_ap09Fluxu_ap09FluxErrumagrmaggmagumagErr
degdegnJynJy
objectint64float32int64float64float64boolint64int64float64float64float64float64float64float64
8769499238_114851411419626224182.9079148514114196262241857.7753558-39.0615181True363445250.601761574.494300925.40253970160046724.0099699650672524.0818583408232526.71969237793988
9706148291_114862846340555160572.83135148628463405551605762.2658894-39.5813315True193448162.203224464.520724725.87485129200129224.3326117849252324.4574696594755126.87575190809596
9762571970_114865661090322262621.98357148656610903222626265.0360865-40.0380382True23449193.337962367.186777725.6842071574959424.3787925708263924.3842910899945126.831790468293185
8769504807_115668040695808055002.91797156680406958080550057.7628398-38.3815605True10363564.493964768.317907126.87620231109101724.1544504374442524.2348499398589426.813663616513505
7957723427_115675957179528017042.97819156759571795280170461.8615102-38.6473312True23637209.076791466.672458925.5992354335742224.1448704467929624.18328773534556226.84013381917539

The cell below creates the cutouts for each of the objects in 5 different bands. By inspecting the cutouts and intrinsic redshift, there are a small number of interlopers (at redshifts $z<2.5$ or so) which visibly show some u-band flux.

In [41]:
filt = ["u", "g", "r", "i", "z"]

for i in range(n_plot):
    fig, axs = plt.subplots(1, len(filt), figsize=(25, 25))

    for j, ax in enumerate(fig.axes):
        
        dataId_deep = {'patch': joined_results['patch'][i].item(),
                       'tract': joined_results['tract'][i].item(),
                       'band': filt[j]}
        cutout = make_image_cutout(service, joined_results['coord_ra'][i].item(),
                                   joined_results['coord_dec'][i].item(),
                                   cutout_size=0.001, imtype='deepCoadd',
                                   dataId=dataId_deep,
                                   filename='cutout_' + filt[j] + '_'
                                   + str(i) + '.fits')

        ax.imshow(cutout.image.array, cmap='gray',
                  vmin=0, vmax=0.7, norm='linear', origin='lower')
        ax.set_xticks([])
        ax.set_yticks([])
        ax.text(.1, .9, filt[j]+'-band z=' + str(joined_results['ts_redshift'][i].item()),
                color='white')
        
        del dataId_deep, cutout
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

Figure 6: Cutout images in $ugriz$ bands (left to right) for 5 of the LBG candidates (top to bottom).

4. Creating calexp cutouts to inspect changing flux for the diaObject¶

Another useful application of the cutout tool is to visually inspect a large number of images of a source identified with difference image analysis (DIA) in order to confirm its variability.

This example is different in that it requires generating cutouts of calexp files around the time the transient is identified, rather than the deepCoadds as used above.

4.1 Query for the source information needed to plot images before and after the transient event¶

This example will use a known diaObject with ID 1250953961339360185 (one identified and studied in DP0.2 tutorial notebook 07a). The cell below will retrieve the information needed to generate a set of images in order of time, including and after the transient event. The TAP query will identify 8 diaSource IDs measured for this diaObject at the ra/dec, and will also retrieve the ccdVisitId of each time the diaSource was observed in a calexp image.

The last line will sort the query return by date, so that they appear in the correct chronological order.

In [42]:
diaObjectId = 1250953961339360185
query = "SELECT ra, decl, diaObjectId, diaSourceId, "\
        "psFlux, psFluxErr, filterName, midPointTai, "\
        "SNR, ccdVisitId, "\
        "scisql_nanojanskyToAbMag(psFlux) AS psAbMag "\
        "FROM dp02_dc2_catalogs.DiaSource "\
        "WHERE diaObjectId = " + str(diaObjectId)\
        + " AND filterName = 'i' "\
        "ORDER BY midPointTai"
print(query)

SELECT ra, decl, diaObjectId, diaSourceId, psFlux, psFluxErr, filterName, midPointTai, SNR, ccdVisitId, scisql_nanojanskyToAbMag(psFlux) AS psAbMag FROM dp02_dc2_catalogs.DiaSource WHERE diaObjectId = 1250953961339360185 AND filterName = 'i' ORDER BY midPointTai
In [43]:
job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)

Job phase is COMPLETED
In [44]:
diaSrc = job.fetch_result().to_table()

The next cell will retrieve calexp cutouts for each of the dia sources in the query, and plot them. Vmin and vmax set a common scaling so that the image scale is not skewed by a bright neighbor, and emphasizes the lower brightness sources such as the dia source. Each figure also is labeled by the time since the first difference image ($\Delta$t in units MJD).

In [45]:
%matplotlib inline

ra = diaSrc['ra'][0] 
dec = diaSrc['decl'][0]
vmin = -200
vmax = 300    
i=0
for src in diaSrc:
    ccdvisitID = src['ccdVisitId']
    visit = str(ccdvisitID)[:-3]
    detector = str(ccdvisitID)[-3:]
    visit = int(visit)
    detector = int(detector)

    dataId_calexp = {'visit':visit, 'detector':detector}
    test = make_image_cutout(service, ra, dec, cutout_size=0.007,
                             imtype='calexp', dataId=dataId_calexp,
                             filename='cutout_'+str(i)+'.fits')

    fig, ax1 = plt.subplots(1, 1, figsize=(2, 2))

    ax1.imshow(test.image.array, cmap='gray', 
               vmin=vmin, vmax=vmax, origin='lower')
    ax1.set_xticks([])
    ax1.set_yticks([])
    ax1.text(.1, .9,r'$\Delta$t=' + str(round(src['midPointTai']-diaSrc['midPointTai'][0],3)),
                color='white',fontsize=12)
    plt.show()
    i = i+1

    del dataId_calexp, test
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

Figure 7: The eight individual $i$-band cutout images for the diaObject.

4.2 Align the calexp images of the transient¶

In the above plots, it is clear that the calexp image orientations are not aligned with each other, because the calexp images are not necessarily all taken in the same position angle on the sky. The cell below will go through each dia source, identify the calexp file in which it was identified, and create cutouts (as above) except now they will all be aligned with each other. At the same time as creating these new aligned fits files for each cutout, the cell will create png files that will be useful for visualizing the transient behavior (to be used in Section 4.3).

The cutout_size variable is in units of degrees, converted from cutout_size_pix using the LSST Science Camera platescale of 0.2 arcseconds per pixel.

In [46]:
i = 0
cutout_size_pix = 131
cutout_size = cutout_size_pix * 0.2 / 3600.0

for src in diaSrc:
    ccdvisitID = src['ccdVisitId']
    visit = str(ccdvisitID)[:-3]
    detector = str(ccdvisitID)[-3:]
    visit = int(visit)
    detector = int(detector)
    dataId_calexp = {'visit': visit, 'detector': detector}
    if i == 0:
        img_ref = make_image_cutout(service, ra, dec, cutout_size=cutout_size,
                                    imtype='calexp', dataId=dataId_calexp)
        minmax = get_minmax_xy(img_ref, cutout_size_pix)
        im_arr = img_ref.image.array[minmax['minx']: minmax['maxx'],
                                     minmax['miny']: minmax['maxy']]
        fig, ax = plt.subplots(1, 1, figsize=(2, 2))
        plt.imshow(im_arr, origin='lower', cmap='gray', vmin=vmin, vmax=vmax)
        ax.set_xticklabels([])
        ax.set_yticklabels([])
        ax.set_xticks([])
        ax.set_yticks([])
        ax.text(.1, .9, r'$\Delta$t='
                + str(round(src['midPointTai']-diaSrc['midPointTai'][0], 2)),
                color='white', fontsize=12)
        figname = os.path.join(tempdir, 'cutout_' + str(i) + '.png')
        if os.path.isfile(figname):
            os.remove(figname)
        plt.savefig(figname)
        plt.show()
        plt.close()
        del minmax, im_arr
    else:
        img = make_image_cutout(service, ra, dec, cutout_size=cutout_size,
                                imtype='calexp', dataId=dataId_calexp)
        img_warped = warp_img(img_ref, img, img_ref.getWcs(), img.getWcs())
        fig, ax = plt.subplots(1, 1, figsize=(2, 2))
        minmax = get_minmax_xy(img_warped, cutout_size_pix)
        im_arr = img_warped.image.array[minmax['minx']: minmax['maxx'],
                                        minmax['miny']: minmax['maxy']]
        plt.imshow(im_arr, origin='lower', cmap='gray', vmin=vmin, vmax=vmax)
        ax.set_xticklabels([])
        ax.set_yticklabels([])
        ax.set_xticks([])
        ax.set_yticks([])
        ax.text(.1, .9, r'$\Delta$t='
                + str(round(src['midPointTai']-diaSrc['midPointTai'][0], 2)),
                color='white', fontsize=12)
        figname = os.path.join(tempdir, 'cutout_' + str(i) + '.png')
        if os.path.isfile(figname):
            os.remove(figname)
        plt.savefig(figname)
        plt.show()
        plt.close()
        del minmax, im_arr, img_warped, img
    i += 1
del img_ref
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

Figure 8: Similar to Figure 7, but for a smaller cutout size and aligning all of the cutout images.

4.3 Make a gif of the aligned calexp cutouts for the transient¶

The cell below will take all the png files created above in the temp folder, and make a gif movie that demonstrates the transient behavior of the diaobject. The gif is by default stored in the same directory as the images it is generated from.

In [47]:
make_gif(tempdir)
display(dimg(data=open(tempdir+'/animation.gif', 'rb').read()))

<IPython.core.display.Image object>

Figure 9: An animated gif centered on a point source (the transient) that brightens then fades.

5. Remove the temp folder¶

As mentioned in Section 1.2.2, empty the temporary folder.

  1. List the contents of the temp folder.
  2. Remove each file in the temp folder, one at a time.
  3. List again the contents of the temp folder to see it is empty.
  4. Delete the temp folder itself.
In [48]:
os.system('ls ' + tempdir)
for filepath in glob.glob(os.path.join(os.getenv('HOME'), tempdir + '/')+"*"):
    os.remove(filepath)
os.system('ls ' + tempdir)
os.system('rmdir ' + tempdir)

animation.gif
cutout-circle.fits
cutout_0.png
cutout_1.png
cutout_2.png
cutout_3.png
cutout_4.png
cutout_5.png
cutout_6.png
cutout_7.png
Out[48]:
0

6. Exercises for the learner¶

  1. Extend the time baseline for the cutouts of the DIA object prior to the transient detection (hint, use the ForcedSourceOnDiaObject table)
  2. Test the failure modes of the cutout tool: what happens if a cutout is requested around coordinates that do not lie inside the linked image?