Notebook

Displaying images using astronomical framework library (lsst.afw.display)¶

Author(s): Brant Robertson (@brantr)
Maintainer(s): Alex Drlica-Wagner (@kadrlica)
Level: Introductory
Last Verified to Run: 2021-09-11
Verified Stack Release: w_2021_33

Learning Objectives:¶

In this tutorial we will

Show how to access the lsst.afw.display routines.
Use the LSST data Butler to access processed data and inspect it visually.

This tutorial is designed to help users get a brief feel for the lsst.afw.display library that enables the visual inspection of data. The lsst.afw library provides an "Astronomical Framework" (afw) while the lsst.daf.* libraries (see, e.g., daf_base) provides a Data Access Framework (daf). Both libraries are used in this tutorial, with the lsst.daf.persistence library used to access a calibrated exposure (calexp) and the lsst.afw.display library used to show the exposure image on the screen.

This tutorial made use of the LowSurfaceBrightness.ipynb StackClub notebook by Alex Drlica-Wagner. More examples of the use of lsst.afw.display can be found in the Stack .

Logistics¶

This notebook is intended to be run on a LARGE instance of the RSP at lsst-lsp-stable.ncsa.illinois.edu or data.lsst.cloud from a local git clone of the StackClub repo.

Setup¶

In [ ]:

# Site, host, and stack version
! echo $EXTERNAL_INSTANCE_URL
! echo $HOSTNAME
! eups list -s | grep lsst_distrib

Step 0) Import Common Python Libraries¶

The matplotlib, numpy, and astropy libraries are widely used Python libraries for plotting, scientific computing, and astronomical data analysis. We will use these packages in common ways below, including the matplotlib.pyplot plotting sublibrary. We also import the warnings library to prevent some routine warning messages from printing to the screen.

In [ ]:

#allow for matplotlib to create inline plots in our notebook
%matplotlib inline                   
import os                            #imports os library
import numpy as np                   #imports numpy with the alias np
import matplotlib.pyplot as plt      #imports matplotlib.pyplot as plt
import warnings                      #imports the warnings library

Let's go ahead and import from astropy the image stretch limits from the familiar zscale() function.

In [ ]:

from astropy.visualization import ZScaleInterval     #This function allows use to use the `zscale()` rescaling limits function familiar from, e.g., DS9, to adjust the image stretch.
zscale = ZScaleInterval()                            #create an alias to the `ZScaleInterval()` function

And let the kernel know that we're happy not to have some useful warnings printed during this tutorial.

In [ ]:

warnings.simplefilter("ignore", category=FutureWarning)   #prevent some helpful but ancillary warning messages from printing during some LSST DM Release calls
warnings.simplefilter("ignore", category=UserWarning)     #prevent some helpful but ancillary warning messages from printing during some LSST DM Release calls

As a last preparatory task, we set the parameters of matplotlib.pyplot to give us a large default size for an image.

In [ ]:

plt.rcParams['figure.figsize'] = (8.0, 8.0)               #set a large default size for our images

Step 1) Loading the LSST DM Stack¶

First, we load some LSST stack libraries to gain access to the image data and visualization routines we'd like to use.

In [ ]:

import lsst.daf.butler      as dafButler      #load the Butler to access data 
import lsst.afw.display     as afwDisplay     #load lsst.afw.display to gain access to image visualization routines.

Step 2) Importing Data to Visualize¶

To display an image, we must first load some data. These data have been processed with the LSST DM Stack, and are organized in a structure that enables us to access them through the DM Stack Butler. For more information on the Butler, see lsst.daf.butler.

In this tutorial, we provide access to two different data sets:

The first data set contains real imaging data from a reprocessing of HSC data available at NCSA in the data directory /datasets/hsc/repo/rerun/RC/v20_0_0_rc1/DM-25349. We access an image from a specific visit (38938) and ccd (32). This happens to be an HSC z-band exposure.
The second data set contains simulated images from the LSST DESC Data Challenge 2 (DC2). These data are available at NCSA in the data directory /datasets/DC2/DR6/Run2.2i/patched/2021-02-10/rerun/run2.2i-calexp-v1/. We access a single image from a specific visit (512055), raft (R20), and detector (76). This happens to be an i-band exposure.

Once we define a string that contains the data directory, we start the Butler instance using the lsst.daf.persistence library alias dafPersist and its Butler class. The Butler object is initialized with a string containing the data directory we wish to access. Running the cell may take a few moments.

With the Butler instance now generated using our data directory, we can retrieve the desired calibrated exposure by telling the butler which filter, CCD, and visit we wish to view. To do this, we definie dictionary with the required information.

In [ ]:

# Location of the DC2 data repository on this site
URL = os.getenv('EXTERNAL_INSTANCE_URL')
if URL.endswith('data.lsst.cloud'): # IDF
    repo = "s3://butler-us-central1-dp01"
elif URL.endswith('ncsa.illinois.edu'): # NCSA
    repo = "/repo/dc2"
else:
    raise Exception(f"Unrecognized URL: {URL}")

dataset='DC2'
collection='2.2i/runs/DP0.1'
dataId = {'visit': 227982, 'raftName': 'R31', 'detector': 129}

In [ ]:

# Create the butler
butler = dafButler.Butler(repo,collections=collection)

# Retrieve the data using the `butler` instance and its function `get()`
calexp = butler.get('calexp', **dataId)

Step 3.1) Use AFWDisplay to Visualize the Image¶

Now, with a Butler instance defined and a calibrated exposure retrieved, we can use lsst.afw.display to visualize the data. The next task is to let AFWDisplay know that we want it to enroll matplotlib as our default display backend. To do this, we use the setDefaultBackend() function. Remember that we made an alias to lsst.afw.display called afwDisplay, so we'll use that to call setDefaultBackend().

In [ ]:

afwDisplay.setDefaultBackend('matplotlib')  # Use lsst.afw.display with the matplotlib backend

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

First create a matplotlib.pyplot figure using plt.figure() -- this will be familiar to anyone with experience using matplotlib.
Then 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 afw_display.
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.
Finally, we can display the image. Do do this, we provide the mtv() method the image member of our calibrated image retrieved by the butler. We can then use plt.show() to display our figure.

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

In [ ]:

plt.figure()                        #create a matplotlib.pyplot figure
display = afwDisplay.Display()      #get an alias to the lsst.afw.display.Display() method
display.scale('asinh', 'zscale')    #set the image stretch algorithm and range
display.mtv(calexp.image)           #load the image into the display
plt.show()                          #show the corresponding pyplot figure

Often you may want to plot two images side-by-side. This can be done by creating matplotlib subplots. Here we plot the image and its mask side by side.

In [ ]:

fig,ax = plt.subplots(1,2,figsize=(14,7))

plt.sca(ax[0]) # set the first axes as current
display1 = afwDisplay.Display(frame=fig)
display1.scale('linear', 'zscale')
display1.mtv(calexp.image)

plt.sca(ax[1]) # set the second axes as current
display2 = afwDisplay.Display(frame=fig)
display2.mtv(calexp.mask)

plt.tight_layout()

It is also possible to plot the mask on top of the image using the calexp.maskedImage. In this case, the image pixel values are plotted in grayscale with the mask values are overplotted in semi-transparent color. We explort the mask plane a bit more in the next section.

In [ ]:

fig = plt.figure()
display = afwDisplay.Display(frame=fig)
display.scale('linear', 'zscale')
display.mtv(calexp.maskedImage)

Congrats! We've plotted an image using lsst.afw.display!

Step 3.2) 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 calexp tutorial for more details). One other component is the mask plane associated with the image. AFWDisplay provides a nice pre-packaged interface for overplotting the mask associated with an 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).

We'll follow the same steps as above to display the image, but we'll add a few modifications

We explicitly set the transparency of the overplotted mask (0 = transparent, 1 = opaque)
We explicitly set the color of the 'DETECTED' mask plane to 'blue' (i.e. all pixels associated with detected objects).
We pass the full calexp object to mtv instead of just the image plane.

In [ ]:

plt.figure()                                     #create a matplotlib.pyplot figure
display = afwDisplay.Display()               #get an alias to the lsst.afw.display.Display() method
display.scale('asinh', 'zscale')             #set the image stretch algorithm and range
display.setMaskTransparency(0)               #set the transparency of the mask plane (0=opaque; 100=transparent)
display.setMaskPlaneColor('DETECTED','pink') #set the color for a single plane in the mask
display.mtv(calexp)                          #load the image and mask plane into the display
plt.show()                                       #show the corresponding pyplot figure

The display object contains more information about the mask planes that can be accessed

In [ ]:

print("Mask plane bit definitions:\n", display.getMaskPlaneColor())  # Print the colors associated to each plane in the mask
print("\nMask plane methods:\n")
help(display.setMaskPlaneColor)

In [ ]:

# A compact way to print some mask information
for maskName, maskBit in calexp.mask.getMaskPlaneDict().items(): 
    print(f'{maskName:18s} (bit=2**{maskBit:02d}): {display.getMaskPlaneColor(maskName)}')

Note that some mask plane bits have their colors unset (i.e., None). These mask planes have colors assigned on the fly at draw-time by mtv. Unfortunately, there isn't a clear way to find out what color was assigned, and this color may change if additional mask planes are added or removed. The best way to have reproducible mask plane colors is to assign an explicit color to each mask plane with setMaskPlaneColor.

Step 3.3) Use AFWDisplay to Visualize a Coadd Image¶

Now we are going to repeat some of the plotting using a coadd data product. First we grab the coadd image and then display it.

In [ ]:

# First use the Butler to get the coadd
coaddId = {'band':'i', 'skymap': 'DC2', 'tract': 4851, 'patch': 29}
coadd = butler.get('deepCoadd', **coaddId)   

In [ ]:

# Plot the coadd; look at all the mask bits!
plt.figure()                          
display = afwDisplay.Display()    
display.scale('asinh', 'zscale')  
display.mtv(coadd)                
plt.show()                            

Wow, look at that mask plane! We can investigate the mask bits as we did for the calexp.

In [ ]:

# A compact way to print some mask information
print("Mask Plane Dict")
for maskName, maskBit in coadd.mask.getMaskPlaneDict().items(): 
    print(f'{maskName:18s} (bit=2**{maskBit:02d}): {display.getMaskPlaneColor(maskName)}')

It turns out that much of the mask plane is coming from the edges of the sensors that go into building the coadd. To remove this mask bit, we make the SENSOR_EDGE and INEXACT_PSF mask planes transparent. The remaining streaky features in the mask plane coming from the bright objects are coming from the REJECTED bit.

In [ ]:

# Let's make the SENSOR_EDGE mask plane transparent
plt.figure()                                     
display = afwDisplay.Display()               
display.scale('asinh', 'zscale') 
display.setMaskTransparency(100,name='SENSOR_EDGE')
display.setMaskTransparency(100,name='INEXACT_PSF')
#display.setMaskTransparency(100,name='REJECTED')
display.mtv(coadd)
plt.show()                                       

In [ ]:

# We can also pan and zoom
plt.figure()                                     
display = afwDisplay.Display()               
display.scale('asinh', 'zscale') 
display.mtv(coadd.image)
display.zoom(16)
display.pan(5150, 18650)
plt.show()

Step 4) More Information about lsst.afw.display¶

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

In [ ]:

method_list = [func for func in dir(display) if callable(getattr(display, func))]
print(method_list)

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 [ ]:

help(display.scale)

In [ ]:

help(display.mtv)

Further Documentation¶

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