Image Display and Manipulation¶

No description has been provided for this image
Contact authors: Alex Drlica-Wagner, Jeff Carlin
Last verified to run: 2024-10-18
LSST Science Pipelines version: Weekly 2024_37
Container Size: medium
Targeted learning level: beginner

Description: Learn to display and manipulate images using the LSST Science Pipelines.

Skills: Display and manipulate images, explore image mask planes, create cutouts and RGB images.

LSST Data Products: Butler calexp and deepCoadd images, coadd objectTable

Packages: lsst.afw.display, lsst.daf.butler, lsst.geom, lsst.afw.image

Credit: This tutorial is based on the AFW_Display_Demo.ipynb notebook originally written by Brant Robertson and maintained by the LSST Stack Club.

More examples of the use of lsst.afw.display can be found in the LSST Science Pipelines documentation.

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 tutorial is designed to help users get a brief feel for the lsst.afw.display library that enables the visual inspection of image data. The lsst.afw library provides an "Astronomical Framework" (afw) while the lsst.daf.* libraries (see, e.g., daf_base) provide a "Data Access Framework" (daf). Both libraries are used in this tutorial, with the lsst.daf.butler library used to access image data and the lsst.afw.display library used to show the exposure image on the screen using the matplotlib backend.

1.1. Package imports¶

Below, the matplotlib.pyplot sublibrary is imported for plotting, and the astropy.wcs package is imported for dealing with World Coordinate Systems (WCS). The matplotlib and astropy libraries are widely used Python libraries for plotting, scientific computing, and astronomical data analysis. The gc (garbage collector) package is imported to help clear the memory of large plots.

The lsst.afw.display and lsst.afw.image modules are loaded to gain access to the image visualization routines we'd like to use, and the lsst.daf.butler and lsst.rsp.get_tap_service modules are used to access data products. The lsst.geom package is used for dealing with sky coordinates.

In [1]:
import matplotlib.pyplot as plt
from astropy.wcs import WCS
from astropy.visualization import make_lupton_rgb
from astropy import units as u
from astropy.coordinates import SkyCoord
import gc

import lsst.afw.display as afwDisplay
from lsst.afw.image import MultibandExposure
from lsst.daf.butler import Butler
from lsst.rsp import get_tap_service
import lsst.geom as geom

1.2. Define functions and parameters¶

1.2.1. Functions¶

Matplotlib stores the data array associated with an image that is plotted. Since the LSST Charge-Coupled Device (CCD) detector images are large (~4k x 4k pixels), this can eventually lead to a memory overflow, which will cause the notebook kernel to die. To mitigate this issue, we define a function to clean up after we plot them.

In [2]:
def remove_figure(fig):
    """
    Remove a figure to reduce memory footprint.

    Parameters
    ----------
    fig: matplotlib.figure.Figure
        Figure to be removed.

    Returns
    -------
    None
    """
    # get the axes and clear their images
    for ax in fig.get_axes():
        for im in ax.get_images():
            im.remove()
    fig.clf()       # clear the figure
    plt.close(fig)  # close the figure
    gc.collect()    # call the garbage collector

Next we define a function to produce a cutout image from a coadd at the user-provided RA, Dec position.

In [3]:
def cutout_coadd(butler, ra, dec, band='r', datasetType='deepCoadd',
                 skymap=None, cutoutSideLength=51, **kwargs):
    """
    Produce a cutout from a coadd at the given ra, dec position.

    Adapted from DC2 tutorial notebook by Michael Wood-Vasey.

    Parameters
    ----------
    butler: lsst.daf.persistence.Butler
        Helper object providing access to a data repository
    ra: float
        Right ascension of the center of the cutout, in degrees
    dec: float
        Declination of the center of the cutout, in degrees
    band: string
        Filter of the image to load
    datasetType: string ['deepCoadd']
        Which type of coadd to load.  Doesn't support 'calexp'
    skymap: lsst.afw.skyMap.SkyMap [optional]
        Pass in to avoid the Butler read.  Useful if you have lots of them.
    cutoutSideLength: float [optional]
        Size of the cutout region in pixels.

    Returns
    -------
    MaskedImage
    """
    radec = geom.SpherePoint(ra, dec, geom.degrees)
    cutoutSize = geom.ExtentI(cutoutSideLength, cutoutSideLength)

    if skymap is None:
        skymap = butler.get("skyMap")

    # Look up the tract, patch for the RA, Dec
    tractInfo = skymap.findTract(radec)
    patchInfo = tractInfo.findPatch(radec)
    xy = geom.PointI(tractInfo.getWcs().skyToPixel(radec))
    bbox = geom.BoxI(xy - cutoutSize // 2, cutoutSize)
    patch = tractInfo.getSequentialPatchIndex(patchInfo)

    coaddId = {'tract': tractInfo.getId(), 'patch': patch, 'band': band}
    parameters = {'bbox': bbox}

    cutout_image = butler.get(datasetType, parameters=parameters,
                              dataId=coaddId)

    return cutout_image

The following function makes cutouts from calexp images instead of deepCoadd images.

In [4]:
def cutout_calexp(butler, ra, dec, visit, detector, cutoutSideLength=51, **kwargs):
    
    """
    Produce a cutout from a calexp at the given ra, dec position.

    Adapted from cutout_coadd which was adapted from a DC2 tutorial
    notebook by Michael Wood-Vasey.

    Parameters
    ----------
    butler: lsst.daf.persistence.Butler
        Helper object providing access to a data repository
    ra: float
        Right ascension of the center of the cutout, in degrees
    dec: float
        Declination of the center of the cutout, in degrees
    visit: int
        Visit id of the calexp's visit
    detector: int
        Detector for the calexp
    cutoutSideLength: float [optional]
        Size of the cutout region in pixels.

    Returns
    -------
    MaskedImage
    """
    
    dataId = {'visit': visit, 'detector': detector}    
    radec = geom.SpherePoint(ra, dec, geom.degrees)
    cutoutSize = geom.ExtentI(cutoutSideLength, cutoutSideLength)    
    calexp_wcs = butler.get('calexp.wcs', **dataId)
    xy = geom.PointI(calexp_wcs.skyToPixel(radec))
    bbox = geom.BoxI(xy - cutoutSize // 2, cutoutSize)
    parameters = {'bbox': bbox}
    cutout_image = butler.get('calexp', parameters=parameters, **dataId)

    return cutout_image

Below we define a function that we will use to produce a color (RGB) image from a set of input images taken with three different filters. We'll use a tool from Astropy (make_lupton_rgb) that helps create RGB images (see this Astropy guide to creating RGB images for more info), and make a convenience function to combine the steps of making an RGB image. (This function was originally written by Aaron Watkins and Nushkia Chamba for their Stack Club Course project exploring surface photometry in the LSST Science Pipelines.)

In [5]:
def create_rgb(image, bgr="gri", stretch=1, Q=10, scale=None):
    """
    Create an RGB color composite image.

    Parameters
    ----------
    image : `MultibandExposure`
        `MultibandExposure` to display.
    bgr : sequence
        A 3-element sequence of filter names (i.e., keys of the exps dict)
        indicating what band to use for each channel. If `image` only has
        three filters then this parameter is ignored and the filters
        in the image are used.
    stretch: int
        The linear stretch of the image.
    Q: int
        The Asinh softening parameter.
    scale: list of 3 floats, each less than 1. (default: None)
        Re-scales the RGB channels.

    Returns
    -------
    rgb: ndarray
        RGB (integer, 8-bits per channel) colour image as an NxNx3 numpy array.
    """

    # If the image only has 3 bands, reverse the order of the bands
    #   to produce the RGB image
    if len(image) == 3:
        bgr = image.filters

    # Extract the primary image component of each Exposure with the
    #   .image property, and use .array to get a NumPy array view.

    if scale is None:
        r_im = image[bgr[2]].array  # numpy array for the r channel
        g_im = image[bgr[1]].array  # numpy array for the g channel
        b_im = image[bgr[0]].array  # numpy array for the b channel
    else:
        # manually re-scaling the images here
        r_im = image[bgr[2]].array * scale[0]
        g_im = image[bgr[1]].array * scale[1]
        b_im = image[bgr[0]].array * scale[2]

    rgb = make_lupton_rgb(image_r=r_im,
                          image_g=g_im,
                          image_b=b_im,
                          stretch=stretch, Q=Q)
    # "stretch" and "Q" are parameters to stretch and scale the pixel values

    return rgb

1.2.2. Parameters¶

Let afwDisplay know that we want it to use matplotlib as our default display backend by setting the setDefaultBackend() function.

Remember that we made an alias to lsst.afw.display called afwDisplay, so we'll use that to call setDefaultBackend().

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

Set the parameters of matplotlib.pyplot to give us a large default size for an image, and set some other parameters to make the default plot style be color-blind friendly and otherwise make plots look nice.

In [7]:
plt.style.use('tableau-colorblind10')
%matplotlib inline
In [8]:
params = {'axes.labelsize': 28,
          'font.size': 24,
          'legend.fontsize': 14,
          'xtick.major.width': 3,
          'xtick.minor.width': 2,
          'xtick.major.size': 12,
          'xtick.minor.size': 6,
          'xtick.direction': 'in',
          'xtick.top': True,
          'lines.linewidth': 3,
          'axes.linewidth': 3,
          'axes.labelweight': 3,
          'axes.titleweight': 3,
          'ytick.major.width': 3,
          'ytick.minor.width': 2,
          'ytick.major.size': 12,
          'ytick.minor.size': 6,
          'ytick.direction': 'in',
          'ytick.right': True,
          'figure.figsize': [8, 8],
          'figure.facecolor': 'White'
          }
plt.rcParams.update(params)

2. Load the data to visualize¶

To display an image, we must first load some data. The DP0.2 data set contains simulated images from the LSST DESC Data Challenge 2 (DC2) that have been reprocessed by the LSST Project using a more recent version of the LSST Science Pipelines. To access these data, we instantiate a Butler directing it to the dp02 data repository configuration and the 2.2i/runs/DP0.2 collection. For more information on the Butler, see lsst.daf.butler and subsequent tutorials in this series.

In [9]:
butler = Butler('dp02', collections='2.2i/runs/DP0.2')

With a Butler instance generated, we can retrieve the desired calibrated exposure by telling the butler which visit and CCD ("detector") we wish to view. To do this, we define a dictionary with the required information. In this case, we access a single image from a specific visit (192350) and detector (175).

In [10]:
dataId = {'visit': 192350, 'detector': 175, 'band': 'i'}
In [11]:
calexp = butler.get('calexp', **dataId)

The calexp that is returned by the Butler in the previous cell is an ExposureF Python object. Exposures are powerful representations of image data because they contain not only the image data, but also a variance image for uncertainty propagation, a bit mask image, and key-value metadata. In the next section, we will use AFWDisplay to visualize the image and mask associated with this Exposure. More documentation on accessing and visualizing an Exposure be found here.

3. Basic image visualization¶

3.1. Use AFWDisplay to visualize the image¶

We are now set to display the image. To do this, we:

  1. Create a matplotlib.pyplot figure using plt.figure() -- this will be familiar to anyone with experience using matplotlib.
  2. Create an alias to the lsst.afw.display.Display method that will allow us to display the data to the screen. This alias will be called display.
  3. Before showing the data on the screen, we have to decide how to apply an image stretch given the data. The algorithm we'll use is asinh -- familiar from SDSS images -- with a range of values set by zscale. To do this, we use the scale() function provided by lsst.afw.display. See the scale() function definition in the interface.py file of the lsst.afw.display library.
  4. Finally, we can display the image. To do this, we provide the mtv() method the calexp.image member of our calibrated image retrieved by the butler. We can then use plt.show() to display our figure.
  5. After we are done creating and displaying the image, we remove the underlying data from memory.

All these tasks are best done within the same notebook cell.

In [12]:
fig = plt.figure()
display = afwDisplay.Display(frame=fig)
display.scale('asinh', 'zscale')
display.mtv(calexp.image)
plt.show()
remove_figure(fig)
No description has been provided for this image

Figure 1: The calexp image, displayed in grayscale with a scale bar at right. Axes labels are in pixels.

3.2. Use matplotlib imshow to visualize the image¶

To see the image axes in sky coordinates instead of pixel coordinates, a simple option is to use astropy's World Coordinate System (WCS) package, along with matplotlib.pyplot's subplot, imshow, and grid functions. Recall that we imported matplotlib.pyplot as plt already, and that we imported the astropy.wcs.WCS constructor as simply WCS. Find more information about imshow and colormaps (cmap).

To do this, we:

  1. Set the figure's projection to be the WCS of the calexp.
  2. Define the extent in pixel coordinates using the bounding box.
  3. Display the calexp image data array using the gray colormap (cmap) and use approximately the same min and max scale values as above.
  4. Add solid white grid lines.
  5. Label the axes, and show the plot.
  6. Remove the underlying data from memory.
In [13]:
fig = plt.figure()
plt.subplot(projection=WCS(calexp.getWcs().getFitsMetadata()))
calexp_extent = (calexp.getBBox().beginX, calexp.getBBox().endX,
                 calexp.getBBox().beginY, calexp.getBBox().endY)
im = plt.imshow(calexp.image.array, cmap='gray', vmin=-200.0, vmax=400,
                extent=calexp_extent, origin='lower')
plt.grid(color='white', ls='solid')
plt.xlabel('Right Ascension')
plt.ylabel('Declination')
plt.show()
remove_figure(fig)
No description has been provided for this image

Figure 2: The same calexp image as in Figure 1, displayed with a different grayscale that increases the number of visible objects. Axes label and grid are in world coordinates.

3.2.1. Use caution when interpreting imshow gridlines!¶

In Figure 2, the more vertical white grid lines represent RA, and they are labled where they intersect the x-axis. In astronomy, the convention is to display images that are north-up. The rotation angle of the calexp displayed in Figure 2 happens to meet this convention, with east/west in RA along the x-axis, and north/south in Dec along the y-axis.

However, the LSST Science Camera will obtain visits at a range of rotation angles. When visualizing a calexp, it is important to understand that it might be rotated such that east/west is along the y-axis instead. Here is a demonstraton of what that might look like.

Define dataId_2 and calexp_2 for a visit that has a significant rotation angle.

In [14]:
dataId_2 = {'visit': 214463, 'detector': 171}
calexp_2 = butler.get('calexp', **dataId_2)

Print the rotation angle.

In [15]:
print(calexp_2.getInfo().getVisitInfo().boresightRotAngle.asDegrees())

289.875617768335

Print the expected center of the image in RA, in hours-minutes-seconds format.

In [16]:
bbox = calexp_2.getBBox()
wcs = calexp_2.getWcs()
center = wcs.pixelToSky(bbox.centerX, bbox.centerY)
c = SkyCoord(ra=center.getRa().asDegrees()*u.degree,
             dec=center.getDec().asDegrees()*u.degree, frame='icrs')
print("Rotated calexp's image center RA in h:m:s format is:", c.ra.hms)
del bbox, wcs, center, c

Rotated calexp's image center RA in h:m:s format is: hms_tuple(h=4.0, m=8.0, s=8.985702947162082)

The expected image center in RA is thus 4h 8m 9s.

Display the calexp_2 with matplotlib imshow and overlay grid lines.

In [17]:
fig = plt.figure()
ax = plt.subplot(projection=WCS(calexp_2.getWcs().getFitsMetadata()))
calexp_extent = (calexp_2.getBBox().beginX, calexp_2.getBBox().endX,
                 calexp_2.getBBox().beginY, calexp_2.getBBox().endY)
ax.set_xlabel('RA (h:m:s)', fontsize=16)
ax.set_ylabel('Dec (deg)', fontsize=16)
ax.coords['ra'].set_ticks(number=12)
ax.coords['dec'].set_ticks(number=12)
im = plt.imshow(calexp_2.image.array, cmap='gray', vmin=-200.0, vmax=400,
                extent=calexp_extent, origin='lower')
plt.grid(color='white', ls='solid')
No description has been provided for this image

Figure 3: A visualization of rotated calexp_2, with grid lines. At first glance, the RA center of the image appears to be about 4h 8m 40s, which is far off the actual RA center of 4h 8m 9s obtained from the image WCS.

CAUTION: Due to image rotation, the RA grid lines are actually the more horizontal lines in this display, but they are still labeled where they intersect the x-axis. Similarly, the Dec grid lines are the more vertical lines, yet still labeled where they intersect the y-axis. Contrast this with Figure 2, where the RA grid lines are more vertical (as expected for typicaly astronomy image display convention of north-up).

Once the orientation of the image is understood, it becomes clearer that the RA center of the image is, indeed, about 4h 8m 9s.

Delete calexp_2 and the associated figure before moving on.

In [18]:
remove_figure(fig)
del ax, im, calexp_extent
del dataId_2, calexp_2

3.3. Use afwDisplay to visualize the image and mask plane¶

The calexp returned by the butler contains more than just the image pixel values (see the Stack Club calexp tutorial for more details). One other component is the mask associated with the image. A mask is composed of a set of "mask planes," 2D binary bit maps corresponding to pixels that are masked for various reasons (see here for more details).

The afwDisplay package maps each bit in the mask plane to a specific display color. We can view this mapping using the code in the following cell. We can also use the setMaskPlaneColor method to change the colors that afwDisplay uses for each mask plane.

Print the colors associated to each plane in the mask.

In [19]:
print("Mask plane bit definitions:\n", display.getMaskPlaneColor())
print("\nMask plane methods:\n")
help(display.setMaskPlaneColor)

Mask plane bit definitions:
 {'BAD': 'red', 'CR': 'magenta', 'EDGE': 'yellow', 'INTERPOLATED': 'green', 'SATURATED': 'green', 'DETECTED': 'blue', 'DETECTED_NEGATIVE': 'cyan', 'SUSPECT': 'yellow', 'NO_DATA': 'orange', 'INTRP': 'green', 'SAT': 'green'}

Mask plane methods:

Help on method setMaskPlaneColor in module lsst.afw.display.interface:

setMaskPlaneColor(name, color=None) method of lsst.afw.display.interface.Display instance
    Request that mask plane name be displayed as color.
    
    Parameters
    ----------
    name : `str` or `dict`
        Name of mask plane or a dictionary of name -> colorName.
    color : `str`
        The name of the color to use (must be `None` if ``name`` is a
        `dict`).
    
        Colors may be specified as any X11-compliant string (e.g.
        `"orchid"`), or by one of the following constants in
        `lsst.afw.display` : `BLACK`, `WHITE`, `RED`, `BLUE`,
        `GREEN`, `CYAN`, `MAGENTA`, `YELLOW`.
    
        If the color is "ignore" (or `IGNORE`) then that mask plane is not
        displayed.
    
        The advantage of using the symbolic names is that the python
        interpreter can detect typos.

Let's plot the image and mask plane side-by-side using matplotlib subplots.

Use plt.sca(ax[0]) to set the first axis as current, and then plt.sca(ax[1]) to switch to the second axis. Using plt.tight_layout() with multi-axis figures helps to avoid axis overlap or excessive white spaces and results in a nicer-looking plot.

In [20]:
fig, ax = plt.subplots(1, 2, figsize=(14, 7))
plt.sca(ax[0])
display1 = afwDisplay.Display(frame=fig)
display1.scale('linear', 'zscale')
display1.mtv(calexp.image)
plt.sca(ax[1])
display2 = afwDisplay.Display(frame=fig)
display2.mtv(calexp.mask)
plt.tight_layout()
plt.show()
remove_figure(fig)
No description has been provided for this image

Figure 4: At left, the same calexp image as in Figures 1 and 2. At right, the mask plane for the calexp illustrates which pixels have a mask value.

The afwDisplay package also provides a nice interface for plotting the mask on top of the image using the calexp.maskedImage. The mask will also be plotted on top of the image if you pass the calexp itself to mtv (as is done later in this notebook).

In [21]:
fig, ax = plt.subplots()
display = afwDisplay.Display(frame=fig)
display.scale('linear', 'zscale')
display.mtv(calexp.maskedImage)
plt.show()
remove_figure(fig)
No description has been provided for this image

Figure 5: The calexp image with the mask plane overlaid.

To investigate the mask in a bit more detail, follow the same steps as above to display the image, but add a few modifications:

  1. Explicitly set the transparency of the overplotted mask to 10% (0 = opaque, 100 = transparent).
  2. Change the color of the 'DETECTED' mask plane from 'blue' (as above) to 'green' (i.e., all pixels associated with detected objects).
  3. Pass the full calexp object to mtv instead of just the image plane.
In [22]:
fig, ax = plt.subplots()
display = afwDisplay.Display(frame=fig)
display.scale('asinh', 'zscale')
display.setMaskTransparency(10)
display.setMaskPlaneColor('DETECTED', 'green')
display.mtv(calexp)
plt.show()
remove_figure(fig)
No description has been provided for this image

Figure 6: Similar to Figure 4, but with the mask transparency reduced to 10% (mostly opaque) and the color representing 'DETECTED' pixels changed from blue to green.

3.4. More information about lsst.afw.display¶

To find more information about lsst.afw.display, print the method list to see what's available. The next cell will print lsst.afw.display methods to the screen.

In [23]:
method_list = [fun for fun in dir(display) if callable(getattr(display, fun))]
print(method_list)

['Buffering', '_Buffering', '_Display__addMissingMaskPlanes', '__class__', '__del__', '__delattr__', '__dir__', '__enter__', '__eq__', '__exit__', '__format__', '__ge__', '__getattr__', '__getattribute__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', 'centroids', 'close', 'delAllDisplays', 'dot', 'erase', 'flush', 'getActiveCallbackKeys', 'getDefaultBackend', 'getDefaultFrame', 'getDisplay', 'getMaskPlaneColor', 'getMaskTransparency', 'image', 'incrDefaultFrame', 'interact', 'line', 'maskColorGenerator', 'mtv', 'pan', 'scale', 'setCallback', 'setDefaultBackend', 'setDefaultFrame', 'setDefaultImageColormap', 'setDefaultMaskPlaneColor', 'setDefaultMaskTransparency', 'setImageColormap', 'setMaskPlaneColor', 'setMaskTransparency', 'show', 'zoom']

If you'd like to learn more about any given function, please see the lsst.afw.display source code.

You can also read the API documentation about the above functions using the Jupyter notebook help() function:

In [24]:
# help(display.scale)
In [25]:
# help(display.mtv)

Congratulations! We've plotted an image in various ways using lsst.afw.display!

4. Extract cutout images¶

When a closer look at images in the vicinity of a star or galaxy is needed, visualizing a small cutout instead of the full image is more effective.

4.1. A cutout from a calexp¶

In a science use case, the coordinates of the object of interest may be known already, but for this example choose pixels near an extended source in the image above.

In [26]:
wcs = calexp.getWcs()
radec = wcs.pixelToSky(2250, 700)
ra, dec = radec.getRa().asDegrees(), radec.getDec().asDegrees()
print(ra, dec)

53.09004683168595 -33.97759584632161

Use the same visit and detector as were used for the dataId in Section 3.

Notice: In order to make a cutout of a calexp, the RA, Dec, visit, and detector must all be specified (recall that visit uniquely identifies band, so band does not need to be explicitly identified). This is different from making a cutout with a deepCoadd (Section 4.2 below), which only requires RA, Dec, and band. This is because there are multiple calexp images for a given coordinate and band, but only one deepCoadd per band.

In [27]:
visit = 192350
detector = 175

Use the cutout_calexp function defined in Section 1 to extract a 301x301 pixel cutout.

Optional: print the help documentation for the function.

In [28]:
# help(cutout_calexp)
In [29]:
my_cutout_calexp = cutout_calexp(butler, ra, dec, visit, detector, cutoutSideLength=301)

Display the cutout from the calexp.

In [30]:
fig, ax = plt.subplots()
display = afwDisplay.Display(frame=fig)
display.scale('asinh', 'zscale')
display.mtv(my_cutout_calexp.image)
plt.show()
remove_figure(fig)
No description has been provided for this image

Figure 7: A cutout from a calexp image, displayed in grayscale, with scale bar at right and axes in pixels. The three brightest objects are galaxies.

In [31]:
del wcs, radec, ra, dec, visit, detector, my_cutout_calexp

4.2. A cutout from a deepCoadd image¶

Up until now, we have been examining an image from a single CCD detector taken during a single visit. For the rest of this notebook, we will switch to examining coadded images made up of multiple exposures that have been combined. Let's start by taking a look at what a full 4k x 4k pixel coadd "patch" image looks like.

We start by grabbing a deepCoadd image and displaying it:

In [32]:
dataId = {'tract': 4431, 'patch': 17, 'band': 'i'}
datasetType = 'deepCoadd'
coadd = butler.get(datasetType, **dataId)
In [33]:
fig, ax = plt.subplots()
display = afwDisplay.Display(frame=fig)
display.scale('asinh', 'zscale')
display.mtv(coadd.image)
plt.show()
remove_figure(fig)
No description has been provided for this image

Figure 8: A deepCoadd patch image, displayed in grayscale with scale bar at right, showing a rich galaxy cluster in the lower-left corner.

The image above is displaying pixel coordinates (note that the coadd patch is part of a larger coadd image called a "tract", so the pixel values do not start at 0,0), but in general it is more useful to be able to select a region based on RA, Dec coordinates. To do this, we'll use the world coordinate system (WCS) object associated with the image.

Extract the WCS solution, which provides the mapping between XY pixel values and sky coordinates, and print it.

In [34]:
wcs = coadd.getWcs()
print(wcs)

FITS standard SkyWcs:
Sky Origin: (55.6521739130, -31.9834710744)
Pixel Origin: (13999, 13999)
Pixel Scale: 0.2 arcsec/pixel

The cluster seems to be centered at about (X, Y) ~ (12500, 8500). Use the "pixelToSky" method of the WCS to get the sky coordinates.

In [35]:
radec = wcs.pixelToSky(12500, 8500)
ra, dec = radec.getRa().asDegrees(), radec.getDec().asDegrees()
print(ra, dec)

55.7506834813934 -32.28892993702273

Now that we have the RA, Dec coordinates of the cluster, we would like to grab a small cutout of the coadded image at this location. To do this, we've defined a user function cutout_coadd at the beginning of the notebook. This function extracts a cutout from a deep coadd image at a given RA, Dec position and desired image size.

In [36]:
help(cutout_coadd)

Help on function cutout_coadd in module __main__:

cutout_coadd(butler, ra, dec, band='r', datasetType='deepCoadd', skymap=None, cutoutSideLength=51, **kwargs)
    Produce a cutout from a coadd at the given ra, dec position.
    
    Adapted from DC2 tutorial notebook by Michael Wood-Vasey.
    
    Parameters
    ----------
    butler: lsst.daf.persistence.Butler
        Helper object providing access to a data repository
    ra: float
        Right ascension of the center of the cutout, in degrees
    dec: float
        Declination of the center of the cutout, in degrees
    band: string
        Filter of the image to load
    datasetType: string ['deepCoadd']
        Which type of coadd to load.  Doesn't support 'calexp'
    skymap: lsst.afw.skyMap.SkyMap [optional]
        Pass in to avoid the Butler read.  Useful if you have lots of them.
    cutoutSideLength: float [optional]
        Size of the cutout region in pixels.
    
    Returns
    -------
    MaskedImage

Call this function to extract a cutout image that is centered on the galaxy cluster.

In [37]:
cutout_image = cutout_coadd(butler, ra, dec, band='i', datasetType='deepCoadd',
                            cutoutSideLength=501)
print("The size of the cutout in pixels is: ", cutout_image.image.array.shape)

The size of the cutout in pixels is:  (501, 501)
In [38]:
fig, ax = plt.subplots()
display = afwDisplay.Display(frame=fig)
display.scale('asinh', 'zscale')
display.mtv(cutout_image.image)
plt.show()
remove_figure(fig)
No description has been provided for this image

Figure 9: A cutout showing a "zoom in" on Figure 7.

5. Overplotting catalog sources¶

Now that we've been able to examine the coadd image and make a cutout, we may be interested in the corresponding catalog of astronomical sources that was extracted from this image by the LSST Science Pipelines. The TAP service is the recommended way to retrieve catalog data for a notebook, and there are examples on how to do this in previous tutorials.

This section demonstrates how to use TAP to retrieve the catalog data within a polygon defined by the corners of the coadd image.

The following cell extracts the XY values of the corners of the cutout image and converts them to RA, Dec, which are then used as spatial constraints in the query below.

In [39]:
wcs = cutout_image.getWcs()

x0 = float(cutout_image.getX0())
y0 = float(cutout_image.getY0())
width = cutout_image.getWidth()
height = cutout_image.getHeight()

xcorners = [x0, x0+width, x0+width, x0]
ycorners = [y0, y0, y0+width, y0+width]

ra_corners = []
dec_corners = []

for i in range(len(xcorners)):
    radec = wcs.pixelToSky(xcorners[i], ycorners[i])
    ra_corners.append(radec.getRa().asDegrees())
    dec_corners.append(radec.getDec().asDegrees())

Instantiate the TAP service.

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

Create a query statement for Objects within the bounding box corners.

Notice: You may be familiar with the BETWEEN ... AND ... format for ADQL spatial queries, and tempted to use it here to select on RA and Dec ranges. However, the Qserv database does not know how to deal with this construction efficiently, so you are much better off constructing a polygon from the corners of the image as below.

In [41]:
query = "SELECT objectId, coord_ra, coord_dec, x, y, tract, patch " + \
        "FROM dp02_dc2_catalogs.Object " + \
        "WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), " +\
        "POLYGON('ICRS', " + str(ra_corners[0]) + ", " + str(dec_corners[0]) + ", " +\
        str(ra_corners[1]) + ", " + str(dec_corners[1]) + ", " +\
        str(ra_corners[2]) + ", " + str(dec_corners[2]) + ", " +\
        str(ra_corners[3]) + ", " + str(dec_corners[3]) + ")) = 1 AND " +\
        "tract = " + str(dataId['tract']) + " AND patch = " + str(dataId['patch'])
print(query)

SELECT objectId, coord_ra, coord_dec, x, y, tract, patch FROM dp02_dc2_catalogs.Object WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), POLYGON('ICRS', 55.767130114294886, -32.30280456542557, 55.73420103348528, -32.302830133388994, 55.73417607075412, -32.2749976201156, 55.76709513053641, -32.274972079508515)) = 1 AND tract = 4431 AND patch = 17

Run the query as an asynchronous job, then fetch the results.

In [42]:
job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)

Job phase is COMPLETED
In [43]:
results = job.fetch_result()
print(len(results))

214

Convert the results to an Astropy table, then take a look at the table.

In [44]:
tab = results.to_table()
tab
Out[44]:
Table length=214
objectIdcoord_racoord_decxytractpatch
degdegpixpix
int64float64float64float64float64int64int64
190994845447018517255.7393808-32.275177512671.7907618747.6994635443117
190994845447018478855.7613844-32.278822112336.98294798681.7890652443117
190994845447018480855.7659405-32.278635112267.64246278685.0823142443117
190994845447018480955.7657282-32.279387912270.88729598671.5347472443117
190994845447018504055.7589613-32.276362212373.81408788726.1058941443117
190994845447018503955.7596823-32.27618512362.83923838729.2837482443117
190994845447018492555.7639926-32.277725912297.26850788701.4806726443117
190994845447018391255.7634528-32.289860712305.7083128483.0549053443117
190994845447018395555.7504241-32.289471112503.95560928490.2627299443117
.....................
190994845447018384755.7618874-32.30166412329.74414028270.6146066443117
190994845447018385655.7592782-32.302142512369.450748262.0420477443117
190994845447018383555.7564491-32.298975212412.44015968319.0969352443117
190994845447018386255.7500326-32.302152312510.1198988261.9994808443117
190994845447018382655.7521728-32.297885512477.48710118338.7732188443117
190994845447018383055.7573277-32.300609312399.10001568289.6690366443117
190994845447018382855.7566968-32.300114612408.69092478298.5832257443117
190994845447020453755.7584044-32.294796812382.61667798394.2817721443117
190994845447020449655.7641292-32.300246312295.6097788296.0987745443117

Display the image cutout, and overplot the sources on the image. Use display buffering to avoid re-drawing the image after each source is plotted.

In [45]:
fig, ax = plt.subplots()
display = afwDisplay.Display(frame=fig)
display.scale('asinh', 'zscale')
display.mtv(cutout_image.image)
with display.Buffering():
    for i in range(len(tab)):
        display.dot('+', tab[i]['x'], tab[i]['y'], ctype=afwDisplay.RED)
        display.dot('o', tab[i]['x'], tab[i]['y'], size=20, ctype='orange')
plt.show()
remove_figure(fig)
No description has been provided for this image

Figure 10: The same deepCoadd cutout as in Figure 8, but with red crosses and orange circles marking the location of Objects.

Note that since the object catalog is assembled from the multi-band coadd images, not every source will necessarily be detected in our i-band image.

6. Plot an RGB composite image¶

We've found a pretty galaxy cluster, but what if we want to know about the colors of the stars and galaxies in that image? To do that, we can extract images taken with three different filters and then assign those images to the RGB channels of a color image.

We will use the MultiBandExposure object type from the LSST Science Pipelines afw.image package. However, one could also create a similar function that simply accepts three separate images without combining them into a MultiBand object. Notice that the MultibandExposure.fromExposures function needs a list of images and filters.

In [46]:
cutout_image_g = cutout_coadd(butler, ra, dec, band='g',
                              datasetType='deepCoadd', cutoutSideLength=701)
cutout_image_r = cutout_coadd(butler, ra, dec, band='r',
                              datasetType='deepCoadd', cutoutSideLength=701)
cutout_image_i = cutout_coadd(butler, ra, dec, band='i',
                              datasetType='deepCoadd', cutoutSideLength=701)
coadds = [cutout_image_g, cutout_image_r, cutout_image_i]
coadds = MultibandExposure.fromExposures(['g', 'r', 'i'], coadds)

Now we are ready to create the RGB (where, in case you're unfamiliar, the RGB denotes "red," "green," and "blue" channels) image. We will assign the g, r, and i-bands to the B, G, and R channels (i.e., g-band, as the bluest filter, will be in the blue channel, r-band in the green channel, and i-band as red). We will use the function create_rgb that we defined at the top of this notebook.

Finally, we will plot two versions of the image to demonstrate how you can change the relative scaling of each RGB channel to produce an image that better highlights certain features. You can experiment with the "scale = []" keyword, or including different bands besides gri.

In [47]:
fig, ax = plt.subplots(figsize=(20, 20), nrows=1, ncols=2)

rgb_original = create_rgb(coadds.image, bgr=['g', 'r', 'i'], scale=None)
ax[0].imshow(rgb_original, origin='lower')
ax[0].set_title('original', fontsize=30)

ax[1].set_title('re-scaled', fontsize=30)
rgb_scaled = create_rgb(coadds.image, bgr=['g', 'r', 'i'],
                        scale=[0.6, 0.7, 1.0])
ax[1].imshow(rgb_scaled, origin='lower')

ax[0].set_axis_off()
ax[1].set_axis_off()
plt.show()
remove_figure(fig)
No description has been provided for this image

Figure 11: Color images built from slightly larger cutouts that are centered on the same coordinates as in Figures 8 and 9. At left, the default RGB scaling is used; at right, a scaling that gives more weight to blue is used, resulting in a cooler palette.

Congratulations - you have now learned how to display images, create cutouts centered on a given position, and make multi-band color images. Enjoy exploring the DP0 images!

7. Exercises for the learner¶

  1. Create coadd and calexp cutouts for the same RA, Dec, and band, and compare them side-by-side.
  2. Write your own function to create a cutout of a color image.

8. Additional Documentation¶

If you'd like some more information on lsst.afw.display, please have a look at the following websites:

  • Info on image indexing conventions.
  • afw.display Doxygen website
  • afw.display GitHub website
  • Getting Started on Image Display (pipelines.lsst.io)
In [ ]: