Background

Although the DEEP2 Survey Team has provided an outstanding data reduction pipeline for DEIMOS slitmask spectroscopy, no analogous package exists for analyzing imaging data. Nevertheless, by using both IRAF and certain IDL procedures distributed with the DEEP2 pipeline, it is possible to reduce DEIMOS imaging data in a straightforward manner.

The following notes describe in general terms the procedure for reducing a set of dithered DEIMOS images to create an image which could be used for measuring astrometric coordinates of targets to be observed with DEIMOS using slitmasks. The procedure treats each DEIMOS CCD independently, then assembles the resulting images into a single mosaic.

Procedure

  1. Create flats. Use either sky frames or dome flats to create flat images.
    1. In IDL, use deimos_read_chip to extract, bias subtract, and apply gain corrections to the individual CCD images within each DEIMOS mosaic image. See example code.
    2. In IRAF, use the flatcombine task to generate individual flats for each chip. See example code.
    3. In IRAF, generate bad pixel masks which exclude the under-illuminated regions of the detectors as located on the skyflats (search for value=1). See example code.
  2. Process science frames. Use the flats you created in step 1 to flatten your images.
    1. In IDL, use deimos_read_chip to extract, bias subtract, apply gain corrections, and divide the appropriate sky/dome flat into the individual CCD images within each DEIMOS mosaic image. See example code.
    2. In IRAF, update image headers to associate the appropriate bad pixel masks with each image:
      hedit *1.fits bpm bpm1.pl add+ ver- sh+ up+
      hedit *2.fits bpm bpm2.pl add+ ver- sh+ up+
      hedit *3.fits bpm bpm3.pl add+ ver- sh+ up+
      hedit *4.fits bpm bpm4.pl add+ ver- sh+ up+
    3. In IRAF, use your preferred cosmic ray cleaning algorithm to remove cosmic ray hits from images.
    4. Check the quality of the assumed gain corrections for each chip by computing the average sky background levels on all images for chip 1, chip 2, etc. Derive correction factors to equalize the sky values and apply the appropriate correction factors to all images.
    5. Optionally subtract sky background from all images to ease processing.
    6. Determine the spatial pixel offsets required to align all images taken with chip 1, and use the imcombine task in IRAF to shift and co-add all chip 1 frames. Repeat for other chips. See example code.
    7. Determine the spatial pixel offsets required to assemble the 4 independent chip images into a single mosaic. In IRAF, use imcombine to assemble the mosaic. See example code.
Note that the method for assembling the individual chips into a single mosaic fails to account for any geometric distortion in the DEIMOS images. Unless an astrometric standard field is observed to derive these distortions, you will probably not have enough stars in your image to derive these distortion corrections. I found that the distortions were sufficiently small that this simple approach was adequate for producing a mosaic image which could then be used to derive astrometry based on USNO and/or Sloan astrometry.


Sample Code

make_skyflat.pro

; make_skyflat.pro / GDW / 2003-Jan-31
;
; Extract individual chips from the images for easier IRAF digestion.
; Also performs the bias subtraction and gain/non-linearity corrections.
;-----------------------------------------------------------------------

indir = '/h/gwirth/gwirth/team-keck/goods'
outdir = '/h/gwirth/gwirth/team-keck/goods/skyflat'
images = [ 'd1230_0085', $
           'd1230_0086', $
           'd1230_0087', $
           'd1230_0088', $
           'd1230_0089', $
           'd1230_0090', $
           'd1230_0091', $
           'd1230_0092', $
           'd1230_0093', $
           'd1230_0095', $
           'd1230_0096', $
           'd1230_0097', $
           'd1230_0098', $
           'd1230_0099']
n_images = n_elements( images)

for chip=1,4 do begin
    for i=0,n_images-1 do begin
        root = images[i]
        in_image = indir+'/'+root+'.fits.gz'
        out_image = outdir+'/'+root+'_'+stringify(chip)+'.fits'
        if file_test( out_image) then begin
            print, out_image, ': already exists; not creating...'
        endif else begin
            print, out_image, ': creating...'
            data = deimos_read_chip( in_image, chip)
            FITS_WRITE, out_image, data
        endelse
    endfor
endfor
end

make_skyflat.cl

# make_skyflat.cl / GDW / 2003-Jan-31
# 
# Combines individual sky images into a flat.  At each pixel, the 11
# highest and 2 lowest values are rejected, leaving only the 3rd lowest value.
# 
# The images are then normalized to unity within the illuminated region and 
# unilluminated regions are set to unity as well to prevent funky artifacts.
#-----------------------------------------------------------------------

print( "Generating Skyflat1...")
flatcombine( "*1.fits", output="Skyflat1", combine="average", reject="minmax",
	ccdtype="", process=no, subsets=no, delete=no, clobber=no, scale="mode",
	statsec="[500:1500,2500:3200]", nlow=2, nhigh=11, nkeep=1, mclip=yes,
	lsigma=3., hsigma=3., rdnoise="0.", gain="1.", snoise="0.", pclip=-0.5,
	blank=1.)

print( "Generating Skyflat2...")
flatcombine( "*2.fits", output="Skyflat2", combine="average", reject="minmax",
	ccdtype="", process=no, subsets=no, delete=no, clobber=no, scale="mode",
	statsec="[500:1500,2500:3200]", nlow=2, nhigh=11, nkeep=1, mclip=yes,
	lsigma=3., hsigma=3., rdnoise="0.", gain="1.", snoise="0.", pclip=-0.5,
	blank=1.)

print( "Generating Skyflat3...")
flatcombine( "*3.fits", output="Skyflat3", combine="average", reject="minmax",
	ccdtype="", process=no, subsets=no, delete=no, clobber=no, scale="mode",
	statsec="[500:1500,2500:3200]", nlow=2, nhigh=11, nkeep=1, mclip=yes,
	lsigma=3., hsigma=3., rdnoise="0.", gain="1.", snoise="0.", pclip=-0.5,
	blank=1.)

print( "Generating Skyflat4...")
flatcombine( "*4.fits", output="Skyflat4", combine="average", reject="minmax",
	ccdtype="", process=no, subsets=no, delete=no, clobber=no, scale="mode",
	statsec="[500:1500,2500:3200]", nlow=2, nhigh=11, nkeep=1, mclip=yes,
	lsigma=3., hsigma=3., rdnoise="0.", gain="1.", snoise="0.", pclip=-0.5,
	blank=1.)

print( "Normalizing images...")
normalize( "Skyflat?.fits", norm=INDEF, sample_section="[*,2100:3800]",
	lower=INDEF, upper=INDEF, imfd="")

print( "Replacing unilluminated regions...")
imreplace( "Skyflat?.fits", value=1., imaginary=0., lower=INDEF, upper=0.5,
	radius=0.)

make_bpm.cl

# make_bpm.pl / GDW / 2003-Feb-06
# 
# Purpose: generate a bad pixel mask based on under-illuminated pixels
# in the skyflats.  Steps include:
#	- flag unity-values pixels as "bad" (1)
#	- use 3x3 boxcar to remove isolated bad pixels
#	- use 9x9 boxcar to expand bad region
#	- combine expanded mask with original one

string skyflat
string bpm
s1 = "temp1.pl"
s2 = "temp2.pl"
s3 = "temp3.pl"
#-----------------------------------------------------------------------
skyflat = "Skyflat1"
bpm = "bpm1.pl"
imexpr( "a==1 ? 1 : 0", s1, skyflat)
boxcar( s1, s2, 3, 3)		# isolated bad pixels become good
imarith( "1", "-", s2, s3)	# bpm -> gpm
imdelete( s1, ver-)
imdelete( s2, ver-)
boxcar( s3, s1, 9, 9)		# good pixels near bad regions get bad
imarith( "1", "-", s1, s2)	# gpm -> bpm
imexpr( "a==1 || b==1 ? 1 : 0", bpm, skyflat, s2)
imdelete( s1, ver-)
imdelete( s2, ver-)
imdelete( s3, ver-)
#-----------------------------------------------------------------------
skyflat = "Skyflat2"
bpm = "bpm2.pl"
imexpr( "a==1 ? 1 : 0", s1, skyflat)
boxcar( s1, s2, 3, 3)		# isolated bad pixels become good
imarith( "1", "-", s2, s3)	# bpm -> gpm
imdelete( s1, ver-)
imdelete( s2, ver-)
boxcar( s3, s1, 9, 9)		# good pixels near bad regions get bad
imarith( "1", "-", s1, s2)	# gpm -> bpm
imexpr( "a==1 || b==1 ? 1 : 0", bpm, skyflat, s2)
imdelete( s1, ver-)
imdelete( s2, ver-)
imdelete( s3, ver-)
#-----------------------------------------------------------------------
skyflat = "Skyflat3"
bpm = "bpm3.pl"
imexpr( "a==1 ? 1 : 0", s1, skyflat)
boxcar( s1, s2, 3, 3)		# isolated bad pixels become good
imarith( "1", "-", s2, s3)	# bpm -> gpm
imdelete( s1, ver-)
imdelete( s2, ver-)
boxcar( s3, s1, 9, 9)		# good pixels near bad regions get bad
imarith( "1", "-", s1, s2)	# gpm -> bpm
imexpr( "a==1 || b==1 ? 1 : 0", bpm, skyflat, s2)
imdelete( s1, ver-)
imdelete( s2, ver-)
imdelete( s3, ver-)
#-----------------------------------------------------------------------
skyflat = "Skyflat4"
bpm = "bpm4.pl"
imexpr( "a==1 ? 1 : 0", s1, skyflat)
boxcar( s1, s2, 3, 3)		# isolated bad pixels become good
imarith( "1", "-", s2, s3)	# bpm -> gpm
imdelete( s1, ver-)
imdelete( s2, ver-)
boxcar( s3, s1, 9, 9)		# good pixels near bad regions get bad
imarith( "1", "-", s1, s2)	# gpm -> bpm
imexpr( "a==1 || b==1 ? 1 : 0", bpm, skyflat, s2)
imdelete( s1, ver-)
imdelete( s2, ver-)
imdelete( s3, ver-)

make_chips.pro

; make_chips.pro / GDW / 2003-Jan-31
;
; Extract individual chips from the images for easier IRAF digestion.
; Also performs the bias subtraction and gain/non-linearity
; corrections, plus division by the normalized skyflat.
;-----------------------------------------------------------------------

indir = '/h/gwirth/gwirth/team-keck/goods'
outdir = '/h/gwirth/gwirth/team-keck/goods/chips'
images = [ 'd1230_0085',  $
           'd1230_0086', $
           'd1230_0087', $
           'd1230_0088', $
           'd1230_0089', $
           'd1230_0090', $
           'd1230_0091', $
           'd1230_0092', $
           'd1230_0093', $
           'd1230_0095', $
           'd1230_0096', $
           'd1230_0097', $
           'd1230_0098', $
           'd1230_0099']
n_images = n_elements( images)

; loop over CCDs...
for chip=1,4 do begin

    flatfile = indir+'/skyflat/Skyflat'+stringify(chip)+'.fits' 
    fits_read, flatfile, flat
    flat = 1./temporary(flat)

    for i=0,n_images-1 do begin
        root = images[i]
        in_image = indir+'/'+root+'.fits.gz'
        out_image = outdir+'/'+root+'_'+stringify(chip)+'.fits'
        print, out_image, ': creating...'
        data = deimos_read_chip( in_image, chip)
        FITS_WRITE, out_image, data*flat
    endfor
endfor
end

do_imcombine.cl

# do_imcombine.cl / GDW / 2003-Feb-01
# 
# Assemble individual exposures for each chip into a single image
# using imcombine task
#------------------------------------------------------------------------
imdelete( "goods?.fits", ver-)
imdelete( "exp?.pl", ver-)
imdelete( "obpm?.pl", ver-)
imcombine( "@inlist1", "goods1", headers="", bpmasks="obpm1",
	rejmasks="", nrejmasks="", expmasks="exp1", sigmas="",
	logfile="STDOUT", combine="average", reject="none", project=no,
	outtype="real", outlimits="", offsets="offsets_1",
	masktype="goodvalue", maskvalue=1., blank=0., scale="mode",
	zero="none", weight="none", statsec="[60:2000,2200:3800]", expname="", 
	lthreshold=INDEF, hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1,
	mclip=yes, lsigma=3., hsigma=7., rdnoise="3", gain="2", snoise="0.",
	sigscale=0.1, pclip=-0.5, grow=0.)
imcombine( "@inlist2", "goods2", headers="", bpmasks="obpm2",
	rejmasks="", nrejmasks="", expmasks="exp2", sigmas="",
	logfile="STDOUT", combine="average", reject="none", project=no,
	outtype="real", outlimits="", offsets="offsets_2",
	masktype="goodvalue", maskvalue=1., blank=0., scale="mode",
	zero="none", weight="none", statsec="[60:2000,2200:3800]", expname="", 
	lthreshold=INDEF, hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1,
	mclip=yes, lsigma=3., hsigma=7., rdnoise="3", gain="2", snoise="0.",
	sigscale=0.1, pclip=-0.5, grow=0.)
imcombine( "@inlist3", "goods3", headers="", bpmasks="obpm3",
	rejmasks="", nrejmasks="", expmasks="exp3", sigmas="",
	logfile="STDOUT", combine="average", reject="none", project=no,
	outtype="real", outlimits="", offsets="offsets_3",
	masktype="goodvalue", maskvalue=1., blank=0., scale="mode",
	zero="none", weight="none", statsec="[60:2000,2200:3800]", expname="", 
	lthreshold=INDEF, hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1,
	mclip=yes, lsigma=3., hsigma=7., rdnoise="3", gain="2", snoise="0.",
	sigscale=0.1, pclip=-0.5, grow=0.)
imcombine( "@inlist4", "goods4", headers="", bpmasks="obpm4",
	rejmasks="", nrejmasks="", expmasks="exp4", sigmas="",
	logfile="STDOUT", combine="average", reject="none", project=no,
	outtype="real", outlimits="", offsets="offsets_4",
	masktype="goodvalue", maskvalue=1., blank=0., scale="mode",
	zero="none", weight="none", statsec="[60:2000,2200:3800]", expname="", 
	lthreshold=INDEF, hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1,
	mclip=yes, lsigma=3., hsigma=7., rdnoise="3", gain="2", snoise="0.",
	sigscale=0.1, pclip=-0.5, grow=0.)

do_imcombine2.cl

# do_imcombine2.cl / GDW / 2003-Feb-01
# 
# Assemble 4 chip images into single DEIMOS mosaic image using
# imcombine task
#------------------------------------------------------------------------
imdel( "goods_full.fits", ver-)
imdel( "obpm_full.pl", ver-)
imdel( "exp_full.pl", ver-)
imcombine( "goods?.fits", 
	"goods_full", headers="", bpmasks="obpm_full",
	rejmasks="", nrejmasks="", expmasks="exp_full", sigmas="",
	logfile="STDOUT", combine="average", reject="none", project=no,
	outtype="real", outlimits="", offsets="offsets_full",
	masktype="badvalue", maskvalue=1., blank=0., scale="mode",
	zero="none", weight="none", statsec="[*,2200:2600]", expname="", 
	lthreshold=INDEF, hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1,
	mclip=yes, lsigma=3., hsigma=7., rdnoise="3", gain="2", snoise="0.",
	sigscale=0.1, pclip=-0.5, grow=0.)