Getting started tutorial part 5: coadding images

In this part of the tutorial series you will combine the individual exposures produced by the singleFrame pipeline (from part 2) into deeper coadds (mosaic images). The dataset that defines how images are reprojected for coaddition is called a skymap. The example repository has a skymap that we can use for this purpose. We will warp (reproject) images into that skymap. Then, you will coadd the warped images together into deep images.

Set up

Pick up your shell session where you left off in part 4. For convenience, start in the top directory of the example git repository.

cd $RC2_SUBSET_DIR

The lsst_distrib package also needs to be set up in your shell environment. See Setting up installed LSST Science Pipelines for details on doing this.

About skymaps

Before you get started, let’s talk about skymaps.

A skymap is a tiling of the celestial sphere, and is used as coordinate system for the final coadded image. A skymap is composed of one or more tracts. Those tracts contain smaller regions called patches. Both tracts and patches overlap their neighbors.

Each tract has a different world coordinate system (WCS), but the WCSs of the patches within a given tract are just linearly-offset versions of the same WCS.

There are two general categories of skymaps:

  1. Whole sky.
  2. A skymap based on the bounding box of a set of input exposures.

Though the HSC dataset you are working with is small, the full dataset is large so we include the skymap computed from the full set and it is of the first type. Using a skymap that describes the full sky has the benefit that you can compare directly with data products from larger data processing runs since the tracts and patches will align exactly.

Warping images onto the skymap

Before assembling the coadded image, you need to warp the exposures created by the singleFrame pipeline onto the pixel grids of patches described in the skymap. You can use the makeWarp pipeline for this.

The data to process can be specified with a dataset query passed to the -d argument, short for --data-query. The example command below only creates coadds for a subset of the data. The remaining sections only use a few patches that have coverage from all the input visits, so downsampling in the coadd generation phase can save time. All the data can be coadded by simply leaving off the -d switch and the arguments to it. Remember that when specifying the output tract and patch information, you must also specify a valid skymap. In this case, the skymap, called hsc_ring_v1 from a larger HSC run is used. The filtering we will do here is:

-d "tract = 9813 AND skymap = 'hsc_rings_v1' AND patch in (38, 39, 40, 41)"

The above dataset query has been tailored to work with the example dataset described in the first tutorial.

The patch in (38, 39, 40, 41) clause selects patches 38, 39, 40 and 41 in the skymap. All five bands for each of the patches will be processed. You can retrieve the skymap and interrogate it using the various member functions.

from lsst.daf.butler import Butler
butler = Butler('SMALL_HSC')
skymap = butler.get('skyMap', skymap='hsc_rings_v1', collections='HSC/RC2/defaults')
tractInfo = skymap.generateTract(9813)
patch = tractInfo[41]
patch.getIndex()
--> (5, 4)

The data queries will typically use the integer id for the patches, but you can use code like that above to find out the x and y indices into the tract’s rectilinear grid of patches. In this case, the patch with id 41 is located at the position (5, 4) in tract 9813.

To warp the images, you will use the pipetask run command again. This time you will specify the makeWarp subset and an appropriate output collection. This example uses coadds as the output collection.

pipetask run -b $RC2_SUBSET_DIR/SMALL_HSC/butler.yaml \
             -d "tract = 9813 AND skymap = 'hsc_rings_v1' AND patch in (38, 39, 40, 41)" \
             -p $RC2_SUBSET_DIR/pipelines/DRP.yaml#makeWarp \
             -i u/$USER/source_calibration,u/$USER/jointcal,u/$USER/fgcm \
             -o u/$USER/warps \
             --register-dataset-types

Note that warping requires the ouptuts of both jointcal and FGCM, so both of those collections need to be specified as inputs. Again, this will warp all calibrated exposures. If you wish to pare down the data to be processed, you can specify a data query like the one earlier in this section using the -d argument.

Tip

As with the singleFrame pipeline, warping only needs the data from an input visit and the skymap. Each warp can be done independently of every other warp. That means it is a good candidate for running in parallel. If you have access to more than one core for processing, specifying the -j argument will speed up this step.

Coadding warped images

Now you’ll assemble the warped images into coadditions for each patch with the assembleCoadd pipeline. As before, we will run without a data query to process a subset of the data, but a selection can be made with the -d argument just as with warping. In this case the the -d arguement could be omitted since the coaddition process will only find the warps from the previous command and will thus only produce coadds for those patches.

Run:

pipetask run -b $RC2_SUBSET_DIR/SMALL_HSC/butler.yaml \
             -d "tract = 9813 AND skymap = 'hsc_rings_v1' AND patch in (38, 39, 40, 41)" \
             -p $RC2_SUBSET_DIR/pipelines/DRP.yaml#assembleCoadd \
             -i u/$USER/warps \
             -o u/$USER/coadds \
             --register-dataset-types

Tip

While coaddition can be done in parallel, each process is more memory intensive than warping because multiple visits from multiple detectors may be put in memory at once. Still, if you have access to a machine with a fair amount of memory, the -j option may still speed up this step.

Wrap up

In this tutorial, you’ve warped exposures into a pre-existing skymap, and then coadded the exposures to make deep mosaics. Here are some key takeaways:

  • Skymaps define the WCS of coadditions.
  • Skymaps are composed of tracts, each of which is composed of smaller patches.
  • The makeWarp pipeline warps exposures into the WCSs of the skymap.
  • The assembleCoadd pipeline coadds warped exposures into deep mosaics.

Continue this tutorial in part 6, where you’ll measure sources in the coadds.