Basic Image Processing Cookbook

General Procedure

In this example, we will refer to a fictitious night of data which, in order of acquisition, looks like this. If you have very short exposure times or very narrow filters (so the sky level is near zero), you may want to follow a more tedious but exact processing procedure.

 

1) If you didn't do this before: from the CCD PC, ftp all frames to your account on baade.

2) Convert fits to imh file formats and put Darktime, ST, Airmass, HJD into headers.

3) Prepare the Zero frame.

4) Create a Dark frame for the correct CCD temperature.

5) Prepare the Flat Field frames for each filter used. Note: changes to come, dome & sky flats.

6) Process all star observations (do Trim, Zero, Dark, and Flat corrections).

7) Notify team members where images can be found.

8) ACLa will backup all processed and raw images to tape, and inform you when images can be deleted.

 

At any point, if you have trouble, you may get a jump-start from the troubleshooting page.

Updated 2007 April 12 - ACL

a FYI, ACL = AndyLayden (a.k.a. YFL, Your Fearless Leader)


FTP Data Transfer

Overview: Copy the images you took from the CCD PC in the warmroom to baade. It is handy to do this immediately after you observe, while the CCD is warming (it only takes a few minutes).

In this example, we will process the data in a directory on baade called /data/layden/PHOT/02MAY30. You should store the unprocessed *.fit images in a subdirectory called RAW (in this example, /data/layden/PHOT/02MAY30/RAW), so if something goes wrong with processing the images you can easily recover the originals. It is also easy for use to archive (copy to tape or DVD) the raw images.

1) First, enter MaxIM_DL (shortcut on desktop) and look through the list of images you took. Do any have incorrect titles (according to our naming convention)? If so, correct them now.

2) From the CCD PC, get into FTP ("file transfer protocol"):

Updated 2006 Aug 08 - ACL


Populate Image Header with Required Keywords

Overview: The image header contains many "keywords" and their values that are required for image processing and analysis (type "imhead yourfile.imh l+" to look at the keywords and values for a file called yourfile.imh). In particular, CCDPROC needs the keyword "DARKTIME" to do the dark current corrections, but MaxIM_DL does not write DARKTIME into the header.

Also, we will express the time of our variable star observations in "Heliocentric Julian Date" (HJD), the number of days since some semi-arbitrary date about 2.5 million days ago. IRAF needs to know the Universal Time (UT), UT-date (DATE-OBS), exposure time (EXPTIME), and Sidereal Time (ST) of the observation in order to compute HJD. We need to estimate the ST when each image was taken from the UT (written into the header) and from the UT/ST "zeropoint" you logged at the beginning of the night.

Note: This and the remaining steps in this manual will be done while logged into baade (remotely).

1) First, there are a bunch of files and programs that you will need for processing. Here, we will copy them all to the current directory on baade. [Hint: Cut/paste this command with the mouse to ensure accuracy.]

cl> cp /data/layden/BGTEL/PROCESS/* .

2) The following IRAF script will convert all of the raw, FITS-format images (RAW/*.fit) into IMH-format images in the current directory (*.imh). [To see the things it has done, type "page readwrite.cl"].

cl> cl < fits2imh.cl
cl> imhead *.imh
(look at all those files!!)

Keep the FITS-format images in /RAW! They are a convenient backup in case you make an error in later processing steps. Plus, we will make a tape copy of the raw and processed images for our archive. ACL will inform you when it is OK to delete the images in /RAW .

3) We will correct all the keywords using a task written by Nick Pearson (1999). This script will ask you several questions, and you should respond as follows (in red).

cl> cl < headadd.cl 
Name of file containing stars RA and Dec: stars.allist  
Epoch of coordinates: 2000
UT at beginning of night: 12 03 14.3 (for example)  
ST at beginning of night: 07 14 00.0 (for example)  

The UT and LST entries are the times you wrote on the top of the logsheet -- they enable the program to compute the ST at any time during the night. Errors made here secretly propagate throughout the analysis procedure and ruin our final light curves, so take care! If you don't have the UT/ST zero-point from the logsheet, look here.

This corrects all the header information. A log of all the operations performed is written to the file "headadd.log", and the file "headadd_go.cl" is the script that was run.

Just to be sure all went well, compare the old values (from MaxIM_DL) with the new values (that IRAF wants). They should be the same (no darktime for biases):

cl> hselect *.imh $I,exptime,darktime yes

Also, look to see that all the star images got reasonable values of UTdate/time, ST, airmass, and HJD (compare with the log sheets you wrote while observing).

cl> hselect *_*.imh $I,date-obs,time-obs,st,airmass,hjd yes

Note: If there are any errors, contact ACL. The program is pretty tricky, and it is not 100% "bulletproof".

Updated 2006 Aug 08 - ACL


Prepare Zero Frame

Overview: The 8 or so bias frames taken at the beginning and/or end of the night will be trimmed (the rows Y>1019 and the columns X>1022 are bad), and then combined into a master frame, "Zero.imh", which contains a 2-D map of the intensity-structure in the bias. The median value of Zero.imh should be about 505 (+/- 5 or so) counts. 

1) Make a list of the bias frames taken at the beginning and/or end of the night (when the dome was shut and the lights were off). These are probably the only biases you have. Check to be sure using IMHEAD. If there are other biases in your list, delete the filenames using emacs [Hint: these emacs keystrokes may be helpful in your editing here and elsewhere].

cl> files *bias*.imh >> bias.list
cl> imhead @bias.list
cl> emacs bias.list (if needed)

2) Open the following packages in IRAF. The file "dpar.ccdred" is a file (page it if you like) that contains the parameter settings for the IRAF task called CCDRED. The "cl < ..." command enters those parameters into the task -- easier than doing it manually with 'epar'. The link shows the correct parameters for CCDRED.

cl> noao
no> imred
im> ccdred
im> cl < dpar.ccdred2
cc> lpar ccdred

3) Use CCDPROC to trim the bad columns off the bias images, and convert them to real format. Use the "cl < ..." command to inject the parameters directly into CCDPROC. Use LPAR to see what the settings are, and run it with the CCDPROC command..

cc> cl < dpar.ccdproc_bias
cc> lpar ccdproc
cc> ccdproc @bias.list

4) Check out the results using CCDLIST. The images are now smaller in size, [1022,1019] rather than [1024,1024]. They are also real (decimal) format rather than short (integer), so they take up more space on the disk. The [T] means that only the Trim processing step was performed on these images.

cc> ccdlist @bias.list

5) Combine these trimmed biases into a high signal-to-noise map of the 2-D structure produced when the chip is read out. The final zero frame is called "Zero.imh". Again, the "cl < ..." command injects our desired parameters directly into ZEROCOMBINE without the need for epar.

cc> cl < dpar.zerocomb
cc> lpar zerocombine
cc> zerocombine

6) Display "Zero.imh" to be sure it looks reasonable. Our zeros are usually pretty uniform everywhere. The first column (X=1) is often significantly higher, and the next few columns (2<X<5) are often a tad higher, but this is normal. The median value of Zero.imh should be around at about 505 +/- 5 counts, and a typical zero image looks like this. Discuss your Zero.imh with ACL if there are any questions.

* If you don't have one already, open an ximtool (image display) window by typing in your IRAF window:
cc> !ximtool &
cc> set stdimage=imt1024
cc> display Zero.imh 1
cc> imstatistics Zero.imh fields="image,midpt"

Updated 2002 May 29 - ACL


Create a Dark Frame

Overview: Every time we take an image, whether there is light falling on the chip or not, a small number of counts appears in each pixel due to "dark current" in the chip (thermal motion of the atoms knocks loose some electrons). The amount of dark current is proportional to the exposure time (i.e., it appears at a constant rate, counts/sec). The amount of dark current varies from pixel to pixel, and tends to be larger when the chip is operated at a higher temperature (faster atomic motions in the silicon chip knock off electrons more frequently).

To correct for the dark current, we have created a number of "dark frames" -- high signal-to-noise images recording the counts in each pixel -- each at a different temperature. In this step, we create a dark frame with the temperature matching the temperature setpoint in our observations by using "linear interpolation" beteween the actual dark images immediately cooler and warmer than our CCD temp. Later, when we run CCDPROC on the star and skyflat images, IRAF will scale this dark frame (darktime=500sec) to the darktime of each image and subtract off the scaled dark image.

1) First, determine the average CCD temperature for your star images. from the temperature recorded in several image headers. The first command creates a file called 'ccdtemp' that contains the CCD temperature from the header of each of your star images. Page the file to see what it looks like. The last command uses the CURFIT task to fit a flat line (order=1) to these data, when you hit 'q' to quit, it will print a mess of stuff to your screen, the last line of which is the mean (average) temperature for the night . Write that number on your worksheet in the box labeled TCCD (1):

cc> hselect *_*.imh UT,CCD-TEMP yes > ccdtemp
cc> !page ccdtemp
cc> curfit ccdtemp function=legendre order=1 power+

2) Next, see which temperatures have dark frames available:

cc> imhead /data/layden/BGTEL/DARKS/DarkT*.imh

3) The dark current increases quickly at hotter temperatures, so we must interpolate between the the dark frames with temperatures immediately above and below our CCD operating temperature. Select from the list in Step 2 the temperatures above and below TCCD: call them T2 (above) and T1 (below) -- write them in the boxes on your worksheet labeled (2b) and (2d).

The fractional contribution of each dark frame to our final dark is then:

4) The next step will copy fractions of each dark image to the current directory, and then add them together to get the final dark that is appropriate to our TCCD. Let's say the dark frame with T1is called 'DarkTmNN.imh', and the one with T2 is called 'DarkTmMM.imh' (to help you remember which files you used, write the values for NN and MM in boxes (2a) and (2c) of your worksheet). Here is how these equations are derived from linear interpolation.

Note: for f1 and f2, you should substitute the numbers you computed in Step 3, and for NN and MM, you should substitute the names of the dark frames with temps below and above TCCD, respectively.

cc> epar imarith [set pixt=r, calct=r, hparams=darktime]  
cc> imarith /data/layden/BGTEL/DARKS/DarkTmNN.imh * f2 D1.imh verb+
cc> imarith /data/layden/BGTEL/DARKS/DarkTmMM.imh * f1 D2.imh verb+
cc> imarith D1.imh + D2.imh LDark.imh verb+
cc> unlearn imarith

5) Let's check to see that all worked properly by checking the median counts in your interpolated dark image:

cc> imstat *D*.imh fields="image,midpt" upper=17000

Write the midpt values for LDark in boxes (3) of your worksheet. Let's check graphically whether this is the correct "linear interpolation" between the dark images immediately warmer and cooler than the CCD temperature you observed at. To do this, edit the file containing the darkcounts vs. temp data for the darks (dark_temp.dat) and add a line at the bottom with the CCD temp and dark count rate (midpt) from LDark (put a "90" in the third column to plot these points with an open circle). Then plot the graph using SuperMongo (SM):

cc> emacs dark_temp.dat &
cc> !sm
: macro read sm.dark_temp
:plt
:quit

You should get a plot in the SM window. The solid symbols and the line connecting them are the sequence of dark images taken at different temperatures, and the open circle is your "interpolated" dark image. It SHOULD lie exactly on the line connecting the solid symbols. If you aren't sure, you can zoom in by editing the sm.dark_temp file and changing the X and Y limits below the line marked CHANGE TO ZOOM, then replot using SM. If the open circle does not fall on the line, something is amiss -- see Andy!

Also, the darktime for LDark should be 500 sec:

cc> hsel *D*.imh $I,darktime yes 

Finally, display the dark image -- it should look something like this:

cc> displa LDark.imh 1

Obtaining the correct dark image is a tricky procedure. If you have any doubts, ask Andy. Here is an example that may clarify the procedure.

 Updated 2006 Aug 08 - ACL


Prepare the Flat Field Image

Overview: The purpose of the flat field images is to remove the variations in sensitivity from one pixel to the next. To do this, we obtain several "skyflat" images of the sky (presumably a source of uniform illumination over the chip) in each filter. We combine the images to remove any cosmic rays or stars, which would produce defects on the final flat field (if, on a image of a certain star, the variable or comparison star fell on one of these defects, the photometry for that star would be incorrect!).

To ensure that the stars change position from one skyflat to the next, we moved the telescope a bit between exposures. However, it is easy to forget that in the excitement of taking skyflats, so we will double-check here that there are no "overlap images" (consecutive images where the stars lie on the same XY positions). We need a bare minimum of 3 skyflats in each filter for the star-rejection techniques in IRAF to work; 5 is better, and 10 or more is best, especially if there are some clouds or other non-uniformities in the sky while you were taking flats.

The ideal skyflat has between 6000 and 8000 counts -- this ensures plenty of signal while allowing good dynamic range above (up to 16,383 counts) and below (down to 0 counts) to record the sensitivities of especially "hot" and "cold" pixels. However, it is sometimes not possible to get such well-exposed frames. Skyflats with count levels between 4000-6000, and between 8000-10,000 are OK; skyflats with 2000-4000 counts will even do in a pinch. There is some "art" involved in knowing whether it is best to add in or leave out extra skyflats of less-than-ideal quality, so if you have any questions or doubts, please discuss them with ACL. This is a critical part of the processing procedure!

1) In the current program, you will have obtained skyflats either with the V or I filter. In this example, we will assume the V filter is being used. Make a list the skyflats using the FILES command, and edit it with emacs.

cc> files *flat*.imh >> sflat.list
cc> emacs sflat.list (remove any incorrect file names)  

2) Process the skyflats. Here we will trim them, zero-correct, and dark-correct them using the task CCDPROC. Ensure all is well with CCDLIST (the processed files should have [TZD]).

cc> cl < dpar.ccdproc_flat
cc> lpar ccdproc
cc> ccdproc @sflat.list zero=Zero.imh
cc> ccdlist @sflat.list

3) New Step (added 2004jul23): The CCD has begun to behave non-linearly. That is, doubling the amount of light does not result in twice the counts for a fixed brightness object. We can correct for this as follows (the coefficients in green may change -- updated 2010feb16 -- if you are working on data taken before this date, see Andy!; this step moved here from (6) on 2007apr12):

cc> imreplace @sflat.list val=99999 low=11500 upp=INDEF
cc> imreplace @sflat.list val=-500 low=INDEF upp=-500
cc> irred
ir> irlincor @sflat.list @sflat.list coeff1=1.0875 coeff2=-0.9173 coeff3=2.4086   
ir> bye
cc> imarith @sflat.list + -38.9171 @sflat.list  calctype=real hparam=''
cc> hedit @sflat.list f_lincor "yes_2009jul" add+ ver-

4) Determine the median counts in each skyflat. This will help you decide which flats to keep and which to reject (see Overview). Also review the exposure times (don't use any shorter than 10sec!).

cc> imstat @sflat.list fields="image,midpt" >> sflat.median
cc> hselect @sflat.list $I,exptime yes >> sflat.exptime
cc> !page sflat.median
cc> !page sflat.exptime

5) Check to be sure there are no "overlap images". Display the skyflats in the order they were obtained and blink them. If the stars move so that they don't overlap on sequential images, all should be well. If some of the images do overlap (e.g., you forgot to move the telescope between them), you can only use one (1) of the overlapped images. Pick the one with the counts closest to 7000, and/or with the exposure time closest to 30 sec (shorter exposures may have a little "shutter effect", longer ones tend to have more faint stars). Ask ACL if you have questions. In our example, I did the following:

* If you don't have one already, open an ximtool (image display) window by typing in your IRAF window:
cc> !ximtool &
cc> set stdimage=imt1024
cc> display sflatI03 1
cc> display sflatI04 2
cc> display sflatI05 3
cc> display sflatI06 4
cc> display sflatI07 5
cc> display sflatI13 6
cc> display sflatI14 7

* blink the images (see ACL if you don't know how; ximtool can hold up to 16 images at once).

* In this example, I found that sflatI05 and sflatI06 were overlapped. I05 had median=7245counts at EXPTIME=26sec, while I06 had median=6134counts at EXPTIME=32sec. Tough call -- I used I05. The rest were OK.

6) Delete the names of the images you do not want to use from sflat.list. Reasons for rejection could include (a) star overlap, (b) low/high counts, (c) exposure time too short. See ACL if you have any doubts!

cc> emacs sflat.list


7) Combine the skyflats into a single, high signal-to-noise, starless final flat, "FlatI.imh" (or "FlatV.imh" if you are working with the V filter). New: the two possible rejection methods have been replaced by one more "robust" method whose parameters (nlow and nhigh) will depend slightly on how many good flats you have:

Nflats nlow

nhigh

<7
1
2
7 to 10
1
3
11 to 16
2
4
17 to 22
2
5
>22
3
6

this will use "minmax" rejection to remove the extreme pixel values, and then "average" to combine the remaining vaules. [If you want to learn more about the way the stars and cosmic rays get rejected, type "help flatcombine" and focus on the rejection algorithms.]

cc> cl < dpar.flatcomb
cc> epar flatcombine [set nlow and nhigh according to the table above]
cc> flatcombine @sflat.list out=FlatI.imh [or FlatV.imh]
cc> ccdlist FlatI.imh

8) Display your final flats to be sure they look reasonable. There should be no obvious stars, and the intensities should drop as you move from the center toward the corners. There may be some "bad pixels" scattered here and there that just won't flatten out. They should look something like this. Discuss the flats with ACL if there are any questions.

cc> display FlatI.imh 1

 Updated 2011 Oct 14 - ACL


Process the Star Images

Overview: Here, we will run CCDPROC on the star (object) images to (a) trim them, (b) subtract "Zero.imh" to remove the bias level (about 505 counts) and the 2-D structure produced when the chip is read out, (c) subtract off the appropriate-temperature dark image, and (d) divide by the appropriate flat field image (e.g., "FlatI.imh"). The resulting images should be completely processed and ready for photometry!

1) Using FILES, create a list of all the star images you took (they should all be through either the V or I filter). For the RR Lyrae program, they should all have an underscore ("_") in them. Use CCDRED to check to be sure no spurious images (biases, flats, Zero, Flat*, Dark*) have crept in. If they have, remove their names from the list using EMACS.

cc> files *_*.imh >> star.list
cc> ccdlist @star.list
cc> emacs star.list  (if needed)  

2) Set the parameters for CCDPROC to do the same processing steps as you did on the flats, then run it. Ensure that the processing has gone well ([TZD] for all star images).

cc> cl < dpar.ccdproc_flat
cc> lpar ccdproc
cc> ccdproc @star.list zero=Zero.imh
cc> ccdlist @star.list

3) New Step (added 2004jul23): The CCD has begun to behave non-linearly. That is, doubling the amount of light does not result in twice the counts for a fixed brightness object. We can correct for this as follows (the coefficients in green may change -- updated 2010feb15 -- if you are working on data taken before this date, see Andy!):

cc> imreplace @star.list val=99999 low=11500 upp=INDEF
cc> imreplace @star.list val=-500 low=INDEF upp=-500
cc> irred
ir> irlincor @star.list @star.list coeff1=1.0875 coeff2=-0.9173 coeff3=2.4086   
ir> bye
cc> imarith @star.list + -38.9171 @star.list  calctype=real hparam=''
cc> hedit @star.list f_lincor "yes_2009jul" add+ ver-

2b) Now set the parameters for CCDPROC and run CCDPROC to apply the flat field to the star image (use "flat = FlatV.imh" if you observed with the V filter). Use CCDLIST to ensure that the processing has gone well ([TZDF] for all star images).

cc> cl < dpar.ccdproc_star
cc> lpar ccdproc
cc> ccdproc @star.list flat=FlatI.imh
cc> ccdlist @star.list

4) Display several of the images to be sure the processing looks reasonable. The sky should appear pretty uniformly gray (i.e., constant brightness). On some images, there are long, curving features -- we think these are caused by reflections inside the telescope, the light of bright stars glinting off a corner or shiny surface. For now, we gotta live with them. Discuss the images with ACL if there are any questions.

cc> display AN_Ser1i2 1

5) Create a file called "processing.doc" in which you record any anomalies or concerns that you have about the processing procedure (eg, you might mention if only 5 flats got used... 8 or 10 would have been better) or with the images you have inspected (eg, changes in brightness of the sky with position on the chip, or nasty "ghost images" ). This is a good place to record any comments that were written on the paper log sheet by the observer (eg, "Clear sky at dusk, great flats!", "Full moon rise at UT=2:15", or "Patchy clouds move in around UT=4:30, thicken till force close at UT=5:40").

cc> emacs processing.doc &

 Updated 2006 Nov 27-- ACL


Notifying the other Team Members, Tape Archives

Once you are satisfied with your processing, email ACL <laydena(at)bgsu(dot)edu> with the directory path of your final images, and a list of all the stars which you observed. I will inspect the images, and if all looks good, I will notify the people who are responsible for doing the photometry on the various stars. These people will copy the images into their own accounts for the photometry.

ACL will make DVD copies of the raw and processed images for the archive. If we need to reprocess or re-photometer the data at some time in the future, we can recover the images from the DVDs.

ACL will tell you when the images can be deleted from your account.

 Updated 2000 Feb 23 - ACL