In this tutorial is described a concrete example usage of the analysis pipeline distributed at the HSC Pipeline Binary Distribution site. Before reading this page, install the pipeline first.

  1. Emplace images
  2. Create data for calibration
  3. Analyze CCDs
  4. Stack
  5. Analyze coadd images
  6. Analyze CCDs again with reference to the coadds
  7. Calibrate catalogs
  8. Calibrate images
  9. (Excursus) Example of a script

Emplace images

Crate a data folder

Create a folder at a place you like. All input data and products will be placed in it. I hereafter assume it ~/HSC.

[user@example ~]$ mkdir ~/HSC
[user@example ~]$ # Cast a spell to initiate this folder
[user@example ~]$ echo 'lsst.obs.hsc.HscMapper' > ~/HSC/_mapper

Download input images

Download images for tutorial from the data list in the pipeline distribution site. When the download completes, expand the archive in an arbitrary directory:

[user@example ~]$ tar xvaf rawdata.tar.xz

Move the images

Emplace the downloaded images in ~/HSC:

[user@example ~]$ # Set up hscPipe
[user@example ~]$ setup-hscpipe
[user@example ~]$ cd rawdata/
[user@example ~]$ # Move the images
[user@example ~]$ hscIngestImages.py ~/HSC  */*.fits --mode=move --create

↑ The destination path (~/HSC) must be an absolute path. Otherwise the command may not work.


The images moved to ~/HSC have been registered in ~/HSC/registry.sqlite3. Look into it to check what images there are.

[user@example ~]$ sqlite3 ~/HSC/registry.sqlite3
sqlite> .header on
sqlite> SELECT visit, filter, field, count(visit)
    FROM raw
    GROUP BY visit, filter, field;
sqlite> .q

If you get something like the following list, the emplacement was successful.


↑If you see smaller figures in the rightmost column, you lack some images.

You can notice that there are DOMEFLATs with visit=907002..907010:2 (“first..last:interval”), and DEN_Cs with visit=905496..905504:2.

Create data for calibration

Create flat data

From the dome flats create flat data. (As for the visit numbers of the dome flats, see above)

[user@example ~]$ reduceFlat.py $HOME/HSC \
    --calib=$HOME/HSC/CALIB_tutorial \
    --output=$HOME/HSC/rerun/calib_tutorial \
    --id visit=907002..907010:2 \
    --detrendId calibVersion=calib_tutorial \
    --config isr.doBias=False \
        isr.doDark=False \
        isr.doFringe=False \
    --do-exec --mpiexec "-n 4"

Note carefully: Values of --calib and calibVersion must match the regular expression [A-Za-z0-9_+-]+. It must not, for example, contain an period.

↑ "-n 4" signifies “ by 4 threads.”

↑ Append “--show config” to the command line to see what can be configured with --config. You will also know the default configuration.

After the data are created, register them.

[user@example ~]$ genCalibRegistry.py \
    --root=$HOME/HSC/CALIB_tutorial \
    --camera=HSC \

We here neglect biases, darks and fringes. You may create them in need with similar commands (reduceBias.py, reduceDark.py, reduceFringe.py)

Analyze CCDs

Before proceeding, you must specify astrometry_net_data in order for the astrometry.

[user@example ~]$ setup astrometry_net_data sdss-dr8

Next, invoke the analysis program that processes CCDs. (As for the range of visit, see above)

[user@example ~]$ reduceFrames.py $HOME/HSC \
    --calib=$HOME/HSC/CALIB_tutorial \
    --rerun=tutorial \
    --id visit=905496..905504:2 \
    --config processCcd.isr.doBias=False \
        processCcd.isr.doDark=False \
        processCcd.isr.doFringe=False \
        doFocus=False $(: for hscPipe 3.2.2 or later) \
    --do-exec --mpiexec "-n 4"
List of outputs

↑ The program will create a directory $HOME/HSC/rerun/RERUN: RERUN being the name specified by --rerun=RERUN. (In this case, $HOME/HSC/rerun/tutorial). Output files are written in this directory. These files will act as inputs of other programs that you will use later.

★Look into the outputs★(UNDER CONSTRUCTION)


Create a skymap

Before mosaicking, create a “tract”. We conventionally call tract a region that is covered with a single WCS (or, that is roughly a field of view).

[user@example ~]$ makeDiscreteSkyMap.py $HOME/HSC \
    --calib=$HOME/HSC/CALIB_tutorial \
    --rerun=tutorial \
    --id visit=905496..905504:2
List of outputs

↑ Tract 0 is created by this command.


If, at this point, latex and dvipng are not installed in the system, install them. If, for instance, you are using Redhat:

[user@example ~]$ sudo yum install dvipng texlive-latex

Then mosaic.

[user@example ~]$ mosaic.py $HOME/HSC \
    --calib=$HOME/HSC/CALIB_tutorial \
    --rerun=tutorial \
    --id visit=905496..905504:2 ccd=0..103 tract=0 \
    --config outputDiag=True \
        outputDir=$HOME/HSC/rerun/tutorial/output_mosaic \
List of outputs


A tract is divided into multiple patches (because the whole tract is too large to deal with)。 To stack a tract, you have to specify patches in the tract to stack. Where you can get information of the patches? The information is in (rerun dir)/deepCoadd/skyMap.pickle.

Look into it.

[user@example ~]$ cd $HOME/HSC/rerun/tutorial/deepCoadd
[user@example deepCoadd]$ python
>>> import pickle
>>> tracts = pickle.load(open("skyMap.pickle", "rb"))
>>> tract = tracts[0] # Select tract 0
>>> print tract.getNumPatches()

↑The result of this command will be:

(12, 12)

↑It means there are (12 columns × 12 rows) patches. (See below.)


Once you see it, return to the directory you were.

>>> exit()
[user@example deepCoadd]$ cd ~
[user@example ~]$

We return from digression. If you want to stack patch (0,0), you specify “patch=0,0”. If you want to stack patch (0,0) and (0,1) respectively, you can specify “patch=0,0^0,1”: ^ being the conjunction. When you want all the patches in a tract, you have to specify patch=0,0^0,1^0,2...

You do not want to do that. Instead, do as follows:

[user@example ~]$ # Save the current field separator.
[user@example ~]$ # (The double quotes are required since IFS is usually spaces.)
[user@example ~]$ oldIFS="$IFS"
[user@example ~]$ # The following generates an array (0,0 0,1 0,2 ... 11,11)
[user@example ~]$ patches=({0..11},{0..11})
[user@example ~]$ # Change the field separator to ^
[user@example ~]$ IFS='^'
[user@example ~]$ # Join the array with the field separator
[user@example ~]$ patches="${patches[*]}"
[user@example ~]$ # Restore the field separator
[user@example ~]$ IFS="$oldIFS"
[user@example ~]$ # (Remember you must not fail to restore IFS)

↑It results in the following $patches:

[user@example ~]$ echo $patches

Now, issue the stack command with the use of $patches

[user@example ~]$ stack.py $HOME/HSC \
    --calib=$HOME/HSC/CALIB_tutorial \
    --rerun=tutorial \
    --id filter=HSC-G tract=0 patch=$patches \
    --selectId visit=905496..905504:2 ccd=0..103 \
    --config \
        doBackgroundReference=False \
        makeCoaddTempExp.bgSubtracted=True \
        assembleCoadd.doMatchBackgrounds=False \
    --clobber-config \
    --do-exec --mpiexec "-n 4"
List of outputs

Analyze coadd images

You use forcedPhotCoadd.py to analize coadd images. However, there is currently a bug that the program times out if we cram many patches into it at a time. Instead of specifying patch=$patches, you have to make a for-loop to deal with the patches one by one.

[user@example ~]$ for i in {0..11},{0..11}
    forcedPhotCoadd.py $HOME/HSC \
        --calib=$HOME/HSC/CALIB_tutorial \
        --rerun=tutorial \
        --id filter=HSC-G tract=0 patch=$i \
        --config references.filter="HSC-G"
List of outputs

Analyze CCDs again with reference to the coadds

Again, the program will times out if we input all at a time as “visit=905496..905504:2”. Make a for-loop to process the visits one by one.

[user@example ~]$ # Note it is not 905496..905504:2
[user@example ~]$ for visit in {905496..905504..2}
    forcedPhotCcd.py $HOME/HSC \
        --calib=$HOME/HSC/CALIB_tutorial \
        --rerun=tutorial \
        --id visit=$visit ccd=0..103 tract=0 \
        --config references.filter="HSC-G"
List of outputs

Calibrate catalogs

[user@example ~]$ calibrateCatalog.py $HOME/HSC \
    --calib=$HOME/HSC/CALIB_tutorial \
    --rerun=tutorial \
    --id visit=905496..905504:2 ccd=0..103 tract=0 \
    --config doApplyCalib=False \
    -j 4 # by 4 threads
List of outputs

Calibrate images

[user@example ~]$ calibrateExposure.py $HOME/HSC \
    --calib=$HOME/HSC/CALIB_tutorial \
    --rerun=tutorial \
    --id visit=905496..905504:2 ccd=0..103 tract=0 \
    -j 4 # by 4 threads
List of outputs

That is all.

(Excursus) Example of a script

Suppose you want a large, whole image instead of chopped-up patches. Since hscPipe is a huge aggregate of Python libraries, you can use the libraries in Python scripts.

As an example, I paste below a script to make a large image from patches.

#!/usr/bin/env python
# encoding: utf-8

# A library of AFW (Application FrameWork), related to images
import lsst.afw.image as afwImage
# A library of AFW, related to geometry
import lsst.afw.geom  as afwGeom

import sys

def main():
    Arrange patches to make a large image, and save it in "large.fits"

    make-large-image PATCHES...
        (e.g. make-large-image *,*.fits)

    if len(sys.argv) <= 1:
        print "Input 1 patch or more."
        return 1

    largeImage = makeLargeImageFromPatches(sys.argv[1:])

    print "Saving (Wait a while)..."

def makeLargeImageFromPatches(filenames):
    @param filenames List of file names
    @return afwImage.ExposureF

    # Loading all image would be dangerous by memory limitation.
    # Load metric information only.
    bboxes = []
    for filename in filenames:
        print "Loading a bounding box:", filename, "..."
            # Get the bounding box of this patch
            # in the parent (= the tract)

    # Calculate the bbox that includes all the bboxes.
    #     Return values of getMin/Max method family are inclusive
    #     (range = [min, max]) cf. getBegin/End (used below)
    x0 = min(bbox.getMinX() for bbox in bboxes)
    y0 = min(bbox.getMinY() for bbox in bboxes)
    x1 = max(bbox.getMaxX() for bbox in bboxes)
    y1 = max(bbox.getMaxY() for bbox in bboxes)

    width  = x1 - x0 + 1
    height = y1 - y0 + 1

    # Prepare a large image
    print "Preparing a large image: ({}, {})".format(width, height)
    exposure = afwImage.ExposureF(width, height)

    # Set the position of this image in the tract
    exposure.setXY0(afwGeom.Point2I(x0, y0))

    # Extract the image data from the exposure
    largeImage = exposure.getMaskedImage()

    # Paste patches into largeImage
    for filename in filenames:
        print "Loading a patch:", filename, "..."

        # Load a patch again
        patch = afwImage.ExposureF(filename)
        bbox = patch.getBBox(afwImage.PARENT)

        # The position of this patch within largeImage:
        #     Return values of getBegin/End method family are one past the end
        #     (range = [begin,end) ) cf. getMin/Max (above)
        xp0 = bbox.getBeginX() - x0
        yp0 = bbox.getBeginY() - y0
        xp1 = bbox.getEndX  () - x0
        yp1 = bbox.getEndY  () - y0

        # Paste the patch (As for the notation in [], see numpy's document)
        largeImage[xp0:xp1, yp0:yp1] = patch.getMaskedImage()[:]

        # Blindly trust that all the patches have the same WCS information etc,
        # and set arbitrary one of them in the large image.
        if not exposure.hasWcs() and patch.hasWcs():
            exposure.setCalib   (patch.getCalib   ())
            exposure.setFilter  (patch.getFilter  ())
            exposure.setWcs     (patch.getWcs     ())

    return exposure

if __name__ == "__main__":

↑Save the script in make-large-image.py, and you can use it as follows:

[user@example ~]$ setup-hscpipe
[user@example ~]$ python make-large-image.py \

↑ Execute it in a machine with more than 16GB of memory. Or, Linux may force a double suicide with the script. Viewing the output large image requires that amount of memory, after all.