Using the LSST Stack Multiband Deblender¶

Owner(s): Fred Moolekamp (@fred3m)
Level: Intermediate
Last Verified to Run: 2021-03-26
Verified Stack Release: v21.0.0

This tutorial is designed to illustrate how to execute the multiband deblender (scarlet) in the LSST stack. This includes a brief introduction to LSST stack objects including:

1. Geometry classes from lsst.geom, such as points and boxes.
2. Higher-level astronomical primitives from lsst.afw, such as the Image, Exposure, and Psf classes.
3. Our core algorithmic Task classes, including those for source detection, deblending, and measurement.

This tutorial is based on Jim Bosch's globular cluster tutorial. However, in it's present state scarlet is unable to process the crowded field (most likely) due to poor initial conditions for the sources in the field. So instead we use a less dense region of the image.

Learning Objectives:¶

After working through this tutorial you should be able to:

1. Configure and run the LSST multiband deblender on a test list of objects;
2. Understand its task context in the DRP pipeline.

Logistics¶

This notebook is intended to be runnable on a large instance of the LSP on lsst-lsp-stable.ncsa.illinois.edu from a local git clone of https://github.com/LSSTScienceCollaborations/StackClub.

Set-up¶

Imports¶

We'll start with some standard imports of both LSST and third-party packages.

Data Set¶

This tutorial can be run with two different data sets:

1. Coadded images made from Subaru Hyper Suprime-Cam (HSC) data in the COSMOS field. We've taken a LSST reprocessing of the HSC-SSP UltraDeep COSMOS field (see this page for information on that reprocessing, and this page for the data), and added simulated stars from a scaled SDSS catalog. The result is a very deep image (deeper than the 10-year LSST Deep-Wide-Fast survey, though not as deep as LSST Deep Drilling fields will be) with both a large number of galaxies and region full of stars.
2. Coadded images from the LSST DESC DC2 simulations. We've taken a fairly characeristic region of the Wide-Fast-Deep survey from LSST DESC DC2 and chosen a region around a galaxy cluster to see how the deblending performs.

We'll be retrieving data using the Butler tool, which manages where various datasets are stored on the filesystem (and can in principle manage datasets that aren't even stored as files, though all of these are). You can find more information on the Butler in a set of tutorials Basics/ButlerTutorial.ipynb and Basics/Gen3ButlerTutorial.ipynb.

We start by creating a Butler instance, pointing it at a Data Repository (which here is just a root directory). Datasets managed by a butler are identified by a dictionary Data ID (specifying things like the visit number or sky patch) and a string DatasetType (such as a particular image or catalog). Different DatasetTypes have different keys, while different instances of the same Dataset Type have different values. All of the datasets we use in this tutorial will correspond to the same patch of sky, so they'll have at least the keys in the dictionary in the next cell (they will also have filter, but with different values).

We can now use those to load a set of grizy coadds, which we'll put directly in a dictionary. The result of each Butler.get call is in this case an lsst.afw.image.Exposure object, an image that actually contains three "planes" (the main image, a bit mask, and a variance image) as well as many other objects that describe the image, such as its PSF and WCS. Note that we (confusingly) use Exposures to hold coadd images as well as true single-exposure images, but combine them into a MultibandExposure, which contains an exposure in each band.

The DatasetType here (deepCoadd_calexp or deepCoadd) is a coadd image on which we've already done some additional processing, such as subtracting the background and setting some mask values, and the extra filter argument gets appended to the Data ID.

Making and displaying color composite images¶

We'll start by just looking at the images, as 3-color composites. We'll use astropy to build those as a nice way to demonstrate how to get NumPy arrays from the MultibandImage objects (the images in coadds). (LSST also has code to make 3-color composites using the same algorithm, and in fact the Astropy implementation is based on ours.) We'll use matplotlib to display the images themselves.

We'll use the following function a few times to display color images. It's worth reading through the implementation carefully to see what's going on.

Notice that we can slice MultibandImage objects (as well a MultibandExposure objects) along the filter dimension using the filter names as indices. Like Exposure objects, MultibandExposure objects have image, mask, and variance properties that contain the image, mask plane, and variance of the Exposure respectively. For now we will only worry about the image property, although internal deblending and measurement algorithms make use of all three objects (when available).

Those images are a full "patch", which is our usual unit of processing for coadds - it's about the same size as a single LSST sensor. That's a bit unweildy (just because waiting for processing to happen isn't fun in a tutorial setting), so we'll reload our dict with sub-images centered on the region of interest.

Note that we can load the sub-images directly with the butler, by appending _sub to the DatasetType and passing a bbox argument. Use sampleBBox to select a sub-region of the image (note that we add a small frame around each blend to include more background regions, which are important for the detection algorithm).

Basic Processing¶

Now we'll try the regular LSST processing tasks, with a simpler configuration than we usually use to process coadds, just to avoid being distracted by complexity. This includes

1. Detection (SourceDetectionTask): given an Exposure, find above-threshold regions and peaks within them (Footprints), and create a parent source for each Footprint.
2. Deblending (ScarletDeblendTask): given a MultibandExposure and a catalog of parent sources, create a child source for each peak in every Footprint that contains more than one peak. Each child source is given a HeavyFootprint, which contains both the pixel region that source covers and the fractional pixel values associated with that source. A SourceDeblendTask is also available using the single band SDSS-HSC deblender that takes a single band Exposure).
3. Measurment (SingleFrameMeasurementTask): given an Exposure and a catalog of sources, run a set of "measurement plugins" on each source, using deblended pixel values if it is a child. Notice that measurement is still performed on single band catalogs, since none of the measurement algorithms work for multiband data.

We'll start by importing these, along with the SourceCatalog class we'll use to hold the outputs.

We'll now construct all of these Tasks before actually running any of them. That's because SourceDeblendTask and SingleFrameMeasurementTask are constructed with a Schema object that records what fields they'll produce, and they modify that schema when they're constructed by adding columns to it. When we run the tasks later, they'll need to be given a catalog that includes all of those columns, but we can't add columns to a catalog that already exists.

To recap, the sequence looks like this:

1. Make a (mostly) empty schema.
2. Construct all of the Tasks (in the order you plan to run them), which adds columns to the schema.
3. Make a SourceCatalog object from the complete schema.
4. Pass the same SourceCatalog object to each Task when you run it.

The first step we'll run is detection, which actually returns a new SourceCatalog object rather than working on an existing one.

Instead, it takes a Table object, which is sort of like a factory for records. We won't use it directly after this, and it isn't actually necessary to make a new Table every time you run MultibandDetectionTask (but you can only create one after you're done adding columns to the schema).

Tasks that return anything do so via a lsst.pipe.base.Struct object, which is just a simple collection of named attributes. The only return values we're interested is sources. That's our new SourceCatalog.

Let's take a quick look at what's in that catalog. First off, we can look at its schema:

Note that this includes a lot of columns that were actually added by the deblend or measurement steps; those will all still be blank (0 for integers or flags, NaN for floating-point columns).

In fact, the only columns filled by SourceDetectionTask are the IDs. But it also attaches Footprint objects, which don't appear in the schema. You can retrieve the Footprint by calling getFootprint() on a row:

Footprints have two components:

• a SpanSet, which represents an irregular region on an image via a list of (y, x0, x1) Spans;
• a PeakCatalog, a slightly different kind of catalog whose rows represent peaks within that Footprint.

You can find a Stack Club notebook devoted to Footprints at SourceDetection/Footprints.ipynb.

If we actually look at the footprints in the catalog we see that some have only a single peak, while others have multiple peaks that need to be deblended.

To display only the pixels contained in the footprint (and not other pixels in the bounding box) we create a MultibandFootprint, which is a HeavyFootprint that contains a SpanSet, PeakCatalog, and flux values for all of the pixels in the SpanSet. In this case the flux is the total measured flux in the image, since no deblending has taken place yet.

It's worth noting that while the peaks can have both an integer-valued position and a floating-point position, they're the same right now; SourceDetectionTask currently just finds the pixels that are local minima and doesn't try to find their sub-pixel locations. That's left to the centroider, which is part of the measurement stage.

Before we can get to that point, we need to run the deblender:

ScarletDeblendTask always returns two catalogs, a templateCatalog that contains the model outputs from scarlet and a place holder for afluxCatalog, which will use the scarlet models as weights to redistribute the flux from the input image (in other words they will contain flux-conserved models). For now only the template catalogs are produced.

The deblender itself sets the parent column for each source, which is 0 for objects with no parent, and all of the columns that begin with deblend_ and also adds new rows to the catalog for each child. It does not remove the parent rows it created those child rows from, and this is intentional, because we want to measure both "interpretations" of the blend family: one in which there is only one object (the parent version) and one in which there are several (the children). Before doing any science with the outputs of an LSST catalog, it's important to remove one of those interpretations (typically the parent one). That can be done by looking at the deblend_nChild and parent fields:

1. parent is the ID of the source from which this was deblended, or 0 if the source is itself a parent.
2. deblend_nChild is the number of child sources this source has (so it's 0 for sources that are themselves children or were never blended).

Together, these define two particularly useful filters:

1. deblend_nChild == 0: never-blended object or de-blended child
2. deblend_nChild == 0 and parent == 0: never-blended object

The first is what you'll usually want to use; the second is what to use if you're willing to throw away some objects (possibly many) because you don't trust the deblender.

The last processing step for our purposes is running measurement, which must be done on each catalog, in each band (if we want measurements for all of them):

Due to an unfortunate bug in the deblender task, the resulting catalogs are not contiguous and we need to copy them into new objects to use them appropriately. This step can be avoided in the near future.

Since we care about deblending (for the sake of this tutorial) lets look at the results from one of the blends displayed above (chosen by hand and specified by idx). We use the HeavyFootprints from the catalog sources that have the same parent to build a model for the entire scene, and to compare the results of the flux conserved and scarlet models. In the process we look at both the scarlet and flux conserved models.

Finally we look at the full models and the residuals. As expected, there are no residuals for the flux conserved model since all of the flux in the image (that is within the footprint) is added to one of the sources. In this particular case that works fine, but in instances where one or more sources were not detected this can cause one source to have its flux contaminated with its neighbor.

Exercises¶

1. Select some othersampleBBox regions and run through the code again, from source detection through measurment and blending displays. Don't foget to change the parent index to view only the children of the correct blend.
2. Look at the configuration options on https://github.com/lsst/meas_extensions_scarlet/blob/master/python/lsst/meas/extensions/scarlet/deblend.py and modify some of them to see the changes in results. For example, set config.symmetric=False to turn off symmetry to see the effect of deblending without symmetry.