Image processing guidelines - LSSTDESC/ReprocessingTaskForce GitHub Wiki

The following guidelines have been written for the CFHT/Megacam instrument and the Abell 2261 / MACJ2243.3-0935 galaxy cluster datasets. Adapting them to another instrument or dataset should be straightforward.

The fist step of the image reduction in DM stack is the processCcd command line task. Assuming that the stack is installed according to the Installing the LSST DM stack and the related packages page, here are the steps to set up the directory structure and run processCcd:

  • Note : The following instructions have been updated to be compatible with recent (November 2016) stack versions

Initial directory structure

mkdir -p rawDownload/A2261
cd rawDownload/A2261
wget --content-disposition -i A2261wgetlist    # This will download the compressed fit images
cd ../..

mkdir input
echo 'lsst.obs.cfht.MegacamMapper' > input/_mapper    # Declare an instrument mapper for the DM butler

The butler is the DM component interacting with the data files.

Set up images in a butler compatible directory structure

 setup pipe_tasks
 ingestImages.py input rawDownload/A2261/*.fz --mode link  # The --mode link will create links instead of copying files

Get or create a configuration file for processCcd

The processCcdConfig.py configuration file should be present in the working directory. It contains the list of customized configuration parameters to be used by the task. A sample configuration file is available here

The block related to the Chebychev background estimation is important to get a good background handling while coadding images

run processCcd.py

The following setup should be issued once for all in the current shell:

setup obs_cfht -k -t your_userid
setup astrometry_net_data -k -t your_userid

We are also using the non default psf determiner, psfex which requires to setup the appropriate package:

  setup meas_extensions_psfex

Launching the task is performed as follows:

processCcd.py input --output output --id visit=806305 ccd=12 --configfile processCcdConfig.py

In this example, the script will process the images corresponding to CCD #12 of exposure #806305. The output catalogs will be created in the output directory following a structure compatible with the butler.

In order to process several images in a single command, one can create a file containing the list of images as the following:

cat AB2261_good_r.list
--id visit=806305 ccd=0..35
--id visit=806306 ccd=0..35
--id visit=806307 ccd=0..35
--id visit=806308 ccd=0..35
--id visit=806416 ccd=0..35
--id visit=806417 ccd=0..35
--id visit=806418 ccd=0..35
--id visit=806419 ccd=0..35

Then the command will be:

 processCcd.py input --output output @AB2261_good_r.list --configfile processCcdConfig.py  --clobber-config -j 4

the parameter -j 4 indicates that the task is going to use 4 processors (1 CCD / processor)

Astrometry

Two astrometry solvers are available. By default the LSST stack is using the algorithm originally developed for HSC. An alternative solver based on astrometry.net can be used by adding the following lines in processConfig.py:

 from lsst.meas.astrom.anetAstrometry import ANetAstrometryTask
 root.calibrate.astrometry.retarget(ANetAstrometryTask)
 root.calibrate.astrometry.solver.sipOrder=4     # Degree of the polynomial used for the sip distortion

In order to solve the Astrometry on Abell 2261, one should download the SDSS DR9 Healpix tiles #176 and 503 (see Getting astrometry reference files)

Coaddition

Coadding images requires several steps:

Create a skymap

This step is now greatly simplified by the use of makeDiscreteSkyMap.py which is able to automatically detect the bounding box containing a set of images

a sample configuration file is available here

The command to launch the task is:

makeDiscreteSkyMap.py output --output output/coadd_dir @AB2261_good_r.list --configfile makeDiscreteSkyMapConfig.py

Identify the list of (tract,patch) overlapping the region of interest

From the makeDiscreteSkyMap.py output (see above):

makeDiscreteSkyMap: tract 0 has corners (261.562, 31.406), (259.395, 31.406), (259.372, 33.255), (261.584, 33.255) (RA, Dec deg) and 9 x 9 patches

we need to get the coordinates of the lower left and upper right corners in order to pass them to the following command:

reportPatches.py output/coadd_dir --config raDecRange="259.395, 31.406, 261.584, 33.255" --id tract=0 patch=0,0 filter=r > patches.txt
  • the --id tract=0 patch=0,0 is meaningless but mandatory (this will probably be improved in the future)

Then we need to modify a little bit the patches.txt file:

 sed -e 's/^/--id filter=r /' patches.txt > patches_r.txt

Warp images to adjust them to the sky map patches

We create a file `AB2261_coadd_good_r.list` containing the following:

--selectId filter=r  visit=806305 
--selectId filter=r  visit=806306 
--selectId filter=r  visit=806307
--selectId filter=r  visit=806308 
--selectId filter=r  visit=806416
--selectId filter=r  visit=806417
--selectId filter=r  visit=806418 
--selectId filter=r  visit=806419

and run the makeCoaddTempExp task:

makeCoaddTempExp.py output --output output/coadd_dir --id filter=r @patches_r.txt @AB2261_coadd_good_r.list --configfile makeCoaddTempExpConfig.py -j 4

Sample configuration file

This will create one warped image for each visit/CCD contributing to a each given patch/tract.

One may have to add the option --timeout 9999999 when running in multiprocessor mode (-j) in order to avoid a timeout error.

Assembling the coadded images

assembleCoadd.py output --output output/coadd_dir @patches_r.txt @AB2261_coadd_good_r.list --configfile assembleCoaddConfig.py

assembleCoaddConfig.py

The following image corresponds to tract=0, patch=5,4:

patch 5,4

Comparison of a galaxy found in tract=1, patch=4,4:

to the same in SDSS:

Multiband Processing

Detailed description of Multiband Processing is available here

The successive steps of multiband processing are:

  • detectCoaddSources.py output/coadd_dir --output output/coadd_dir --id filter=u tract=0 patch=8,3
  • mergeCoaddDetections.py output/coadd_dir --output output/coadd_dir --id tract=1 patch=8,3 filter=u^g^r^i^z
  • measureCoaddSources.py output/coadd_dir --output output/coadd_dir --id tract=1 patch=8,3 filter=z
  • mergeCoaddMeasurements.py output/coadd_dir --output output/coadd_dir --id tract=1 patch=8,3 filter=u^g^r^i^z

The configuration files corresponding to each step can be copied from here

The final step is the forcedPhotometry which can be run at the CCD level (forcedPhotCcd.py) or at the coadd level (forcedPhotCoadd.py).

CModel should be turned on in measureCoaddSources.py using the following parameters:

import lsst.meas.modelfit
import lsst.shapelet
config.measurement.plugins.names |= ["modelfit_ShapeletPsfApprox", "modelfit_CModel"]

The galaxy shape measurement using meas_extensions_shapeHSM can be turned on in the multiband processing. In order to do so, the meas_extensions_shapeHSM should be installed first following the special procedure described here

Color diagram of galaxy candidates from patch 2,3 of MACSJ2243.3-0935 where the galaxy fluxes has been determined coherently in the various filters using the CModel fit MACSJ2243.3-0935