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.
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 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
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.
visit|filter|field|count(visit) 905496|HSC-G|DEN_C|112 905498|HSC-G|DEN_C|112 905500|HSC-G|DEN_C|112 905502|HSC-G|DEN_C|112 905504|HSC-G|DEN_C|112 907002|HSC-G|DOMEFLAT|112 907004|HSC-G|DOMEFLAT|112 907006|HSC-G|DOMEFLAT|112 907008|HSC-G|DOMEFLAT|112 907010|HSC-G|DOMEFLAT|112
↑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.
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 \
--create
We here neglect biases, darks and fringes. You may create them in need with similar commands (reduceBias.py, reduceDark.py, reduceFringe.py)
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"
↑ 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)
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
↑ 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 \
--clobber-config
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.)
| (0,0) | (1,0) | ⋯ | (11,0) |
| (0,1) | (1,1) | ⋯ | (11,1) |
| ⋮ | ⋮ | ⋱ | ⋮ |
| (0,11) | (1,11) | ⋯ | (11,11) |
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 0,0^0,1^0,2^0,3^(...omitted...)^11,10^11,11
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"
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}
do
forcedPhotCoadd.py $HOME/HSC \
--calib=$HOME/HSC/CALIB_tutorial \
--rerun=tutorial \
--id filter=HSC-G tract=0 patch=$i \
--config references.filter="HSC-G"
done
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}
do
forcedPhotCcd.py $HOME/HSC \
--calib=$HOME/HSC/CALIB_tutorial \
--rerun=tutorial \
--id visit=$visit ccd=0..103 tract=0 \
--config references.filter="HSC-G"
done
[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
[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
That is all.
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)..." largeImage.writeFits("large.fits") 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, "..." bboxes.append( # Get the bounding box of this patch # in the parent (= the tract) afwImage.ExposureF(filename).getBBox(afwImage.PARENT) ) # 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.setMetadata(patch.getMetadata()) exposure.setWcs (patch.getWcs ()) return exposure if __name__ == "__main__": 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 \
$HOME/HSC/rerun/tutorial/deepCoadd/HSC-G/0/*.fits
↑ 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.