Astronomy‎ > ‎

IRAF

Here are some useful links to check out in preparation for or while processing WIYN mosaic data.
WIYN 0.9 Observing - This site will be familiar if you've observed at Kitt Peak.  However, if you have never been or have no experience with how the data is taken, this might be a helpful site.  It will also get you familar with terms and concepts.
The Mosaic CCD Project - This site is one of the most useful sites.  It has links to telescope manuals, filter information, IRAF processing packages, newsletters, and even some images!

When I first started to learn IRAF and figure out how to process my data, these were the cookbooks/helpguides I used.
A User's Guide to CCD Reductions with IRAF - By Phil Massey.  This is a link to a ps file download.  This was singlehandely the most helpful help guide I used.
The NOAO Deep Wide-Field Survey MOSAIC Data Reductions - Another excellent cookbook.  Matches up excellently with all of my steps.  Also explains the theory behind what is going on much much better than I do.
NOAO MOSAIC-1 CCD IMAGER - This site is not quite a cookbook, but explains theory and some of the electronic equiptment at the telescope.


The following is step-by-step cookbook/notes from my experience reducing WIYN 0.9m mosaic data using IRAF.
Notes also available here to download as a PDF.

last updated: 2009


The following is by no means an exhaustive procedure.  It is meant to be used in conjunction with other sources for a complete understanding.  These steps are also not the right parameters for every situation, merely as a guide based on what I used in my data processing.

 

Tips for beginners:

- It’s always a good idea to examine your images before you start anything.  Make sure they look OK, that you aren’t combining images that look funny, have striping, a chip that has fallen out, etc.

- It’s also a good idea to check the observer’s log.  That will tell you if one of your images is bad as well. 

- Backing up your data at crucial steps is always a life saver.  If you make one mistake at the end of the processing, you won’t have to go back and start all over.

- Create your own log file of what you do to your data, your parameters, problems you have, etc

_________________________________________________________________________________________________

Part 1: Basic Image Flat Fielding/Processing
ccdproc
zerocombine
ccdproc
flatcombine
ccdproc
[skyflat]

 


ccdproc:

First pass of the task ccdproc.  The first pass will be applying cross talk correction, overscan, and trim to all of your images.  In the images parameter you should specify all your types of images.  In my folder I have my dflats, sflats, zeros, and objects.  Otherwise, all of my other parameters were as follows
                PACKAGE = mscred
                   TASK = ccdproc

                images  = dflat*.fits,sflat*.fits,zero*.fits,obj*.fits  List of Mosaic CCD imag
                (output =                             ) List of output processed images
                (bpmasks=                          ) List of output bad pixel masks
                (ccdtype=                            ) CCD image type to process
                (noproc =                      no) List processing steps only?
                (xtalkco=                     yes) Apply crosstalk correction?
                (fixpix =                         no) Apply bad pixel mask correction?
                (oversca=                    yes) Apply overscan strip correction?
                (trim   =                       yes) Trim the image?
                (zerocor=                     no) Apply zero level correction?
                (darkcor=                     no) Apply dark count correction?
                (flatcor=                       no) Apply flat field correction?
                (sflatco=                      no) Apply sky flat field correction?
                (split  =                        no) Use split images during processing?
                (merge  =                    no) Merge amplifiers from same CCD?
                (xtalkfi=              !xtalkfil) Crosstalk file
                (fixfile=                   BPM) List of bad pixel masks
                (saturat=              INDEF) Saturated pixel threshold
                (sgrow  =                       0) Saturated pixel grow radius
                (bleed  =               INDEF) Bleed pixel threshold
                (btrail =                       20) Bleed trail minimum length
                (bgrow  =                      0) Bleed pixel grow radius
                (biassec=           !biassec) Overscan strip image section
                (trimsec=          !trimsec) Trim data section
                (zero   =                   Zero) List of zero level calibration images
                (dark   =                   Dark) List of dark count calibration images
                (flat   =                  Dflat*) List of flat field images
                (sflat  =                  Sflat*) List of secondary flat field images
                (minrepl=                     1.) Minimum flat field value
                (interac=                     no) Fit overscan interactively?
                (functio=             spline3) Fitting function
                (order  =                      15) Number of polynomial terms or spline pieces
                (sample =                     *) Sample points to fit
                (naverag=                     1) Number of sample points to combine
                (niterat=                       1) Number of rejection iterations
                (low_rej=                   3.) Low sigma rejection factor
                (high_re=                   3.) High sigma rejection factor
                (grow   =                    0.) Rejection growing radius
                (fd     =                           )
                (fd2    =                          )
                (mode   =                   ql)


zerocombine:

The next thing to do is combine your zeros together to make one zero file that will be applied to the all the other images.  This task will output you that zero file to apply to others.  Input - your zero file names, output - what your new file will be called.

                PACKAGE = mscred
                  TASK = zerocombine

                input   =                zero*    List of zero level images to combine
                (output =                 Zero) Output zero level name
                (combine=        average) Type of combine operation
                (reject =             crreject) Type of rejection
                (ccdtype=                zero) CCD image type to combine
                (process=                  yes) Process images before combining?
                (delete =                    no) Delete input images after combining?
                (scale  =                 none) Image scaling
                (statsec=                         ) Image section for computing statistics
                (nlow   =                        1) minmax: Number of low pixels to reject
                (nhigh  =                        1) minmax: Number of high pixels to reject
                (nkeep  =                       1) Minimum to keep (pos) or maximum to reject (neg
                (mclip  =                    yes) Use median in sigma clipping algorithms?
                (lsigma =                      3.) Lower sigma clipping factor
                (hsigma =                     3.) Upper sigma clipping factor
                (rdnoise=        RDNOISE) ccdclip: CCD readout noise (electrons)
                (gain   =                 GAIN) ccdclip: CCD gain (electrons/DN)
                (snoise =                     0.) ccdclip: Sensitivity noise (fraction)
                (pclip  =                   -0.5) pclip: Percentile clipping parameter
                (blank  =                     0.) Value if there are no pixels
                (mode   =                   ql)


ccdproc:

This is our second pass with ccdproc.  Using the zero we have just created in zerocombine, we will make the zero level correction on all our other images.  Change your input to all of your images you have left, dflats, sflats, and objects.  Make sure your “zero= Zero” actually matches up with what called your zero image that you just created.
                PACKAGE = mscred
                  TASK = ccdproc

                images  =   dflat*,sflat*,obj*  List of Mosaic CCD images to process
                (output =                      ) List of output processed images
                (bpmasks=                    ) List of output bad pixel masks
                (ccdtype=                     ) CCD image type to process
                (noproc =                no) List processing steps only?
                (xtalkco=                 no) Apply crosstalk correction?
                (fixpix =                   no) Apply bad pixel mask correction?
                (oversca=                no) Apply overscan strip correction?
                (trim   =                     no) Trim the image?
                (zerocor=                yes) Apply zero level correction?
                (darkcor=                 no) Apply dark count correction?
                (flatcor=                   no) Apply flat field correction?
                (sflatco=                   no) Apply sky flat field correction?
                (split  =                     no) Use split images during processing?
                (merge  =                 no) Merge amplifiers from same CCD?
                (xtalkfi=           !xtalkfil) Crosstalk file
                (fixfile=                 BPM) List of bad pixel masks
                (saturat=            INDEF) Saturated pixel threshold
                (sgrow  =                     0) Saturated pixel grow radius
                (bleed  =             INDEF) Bleed pixel threshold
                (btrail =                     20) Bleed trail minimum length
                (bgrow  =                    0) Bleed pixel grow radius
                (biassec=        !biassec) Overscan strip image section
                (trimsec=       !trimsec) Trim data section
                (zero   =                Zero) List of zero level calibration images
                (dark   =                Dark) List of dark count calibration images
                (flat   =               Dflat*) List of flat field images
                (sflat  =               Sflat*) List of secondary flat field images
                (minrepl=                  1.) Minimum flat field value
                (interac=                  no) Fit overscan interactively?
                (functio=          spline3) Fitting function
                (order  =                   15) Number of polynomial terms or spline pieces
                (sample =                   *) Sample points to fit
                (naverag=                   1) Number of sample points to combine
                (niterat=                     1) Number of rejection iterations
                (low_rej=                   3.) Low sigma rejection factor
                (high_re=                   3.) High sigma rejection factor
                (grow   =                    0.) Rejection growing radius
                (fd     =                          )
                (fd2    =                         )
                (mode   =                   ql)


flatcombine:

This task is used to combine flats.  First, we will do our dome flats.  Input your dome flats, and the output will be your combined dome flat image that will be used next.  The task should be smart enough to only combine images of the same band (will create an r-band dome flat from the r-band exposures, a different ha dome flat from the ha exposures).  Remember to check exposure times in the objects’ image headers.  If you have two different exposure times for the same band, you may want to exclude the exposure time that does not match the rest.  You should end up with a Dflat for each band.

                PACKAGE = mscred

                  TASK = flatcombine

                input   =               dflat*  List of flat field images to combine
                (output =               Dflat) Output flat field root name
                (combine=       average) Type of combine operation
                (reject =           avsigclip) Type of rejection
                (ccdtype=                 flat) CCD image type to combine
                (process=                 yes) Process images before combining?
                (subsets=                 yes) Combine images by subset parameter?
                (delete =                    no) Delete input images after combining?
                (scale  =               mode) Image scaling
                (statsec=                        ) Image section for computing statistics
                (nlow   =                      1) minmax: Number of low pixels to reject
                (nhigh  =                      1) minmax: Number of high pixels to reject
                (nkeep  =                     1) Minimum to keep (pos) or maximum to reject (neg
                (mclip  =                  yes) Use median in sigma clipping algorithms?
                (lsigma =                    3.) Lower sigma clipping factor
                (hsigma =                   3.) Upper sigma clipping factor
                (rdnoise=         rdnoise) ccdclip: CCD readout noise (electrons)
                (gain   =                  gain) ccdclip: CCD gain (electrons/DN)
                (snoise =                    0.) ccdclip: Sensitivity noise (fraction)
                (pclip  =                  -0.5) pclip: Percentile clipping parameter
                (blank  =                    1.) Value if there are no pixels
                (mode   =                   q)


ccdproc:

Third pass of ccdproc.  We are now applying the dome flats to our remaining objects and skyflats.  Make sure your “flat= Dflat*” actually matches up with what your dome flats are called.
                PACKAGE = mscred
                   TASK = ccdproc

                images  =          obj*,sflat*  List of Mosaic CCD images to process
                (output =                      ) List of output processed images
                (bpmasks=                    ) List of output bad pixel masks
                (ccdtype=                     ) CCD image type to process
                (noproc =                 no) List processing steps only?
                (xtalkco=                  no) Apply crosstalk correction?
                (fixpix =                    no) Apply bad pixel mask correction?
                (oversca=                 no) Apply overscan strip correction?
                (trim   =                    no) Trim the image?
                (zerocor=                 no) Apply zero level correction?
                (darkcor=                 no) Apply dark count correction?
                (flatcor=                  yes) Apply flat field correction?
                (sflatco=                   no) Apply sky flat field correction?
                (split  =                     no) Use split images during processing?
                (merge  =                 no) Merge amplifiers from same CCD?
                (xtalkfi=           !xtalkfil) Crosstalk file
                (fixfile=                 BPM) List of bad pixel masks
                (saturat=            INDEF) Saturated pixel threshold
                (sgrow  =                     0) Saturated pixel grow radius
                (bleed  =             INDEF) Bleed pixel threshold
                (btrail =                     20) Bleed trail minimum length
                (bgrow  =                    0) Bleed pixel grow radius
                (biassec=        !biassec) Overscan strip image section
                (trimsec=       !trimsec) Trim data section
                (zero   =               Zeros) List of zero level calibration images
                (dark   =                 Dark) List of dark count calibration images
                (flat   =                Dflat*) List of flat field images
                (sflat  =                Sflat*) List of secondary flat field images
                (minrepl=                   1.) Minimum flat field value
                (interac=                   no) Fit overscan interactively?
                (functio=           spline3) Fitting function
                (order  =                    15) Number of polynomial terms or spline pieces
                (sample =                    *) Sample points to fit
                (naverag=                    1) Number of sample points to combine
                (niterat=                      1) Number of rejection iterations
                (low_rej=                   3.) Low sigma rejection factor
                (high_re=                   3.) High sigma rejection factor
                (grow   =                    0.) Rejection growing radius
                (fd     =                           )
                (fd2    =                          )
                (mode   =                   ql)

 

[skyflats]:

In my processing, I excluded using skyflats.  What ever I ended up doing to my images with sky flats made the objects look worse.  If you plan to do sky flats, read more online.  You may have to deal with pupiling problems.  I did not have pupiling to deal with; I combined my sflats with the task sflatcombine, and then ran the task ccdproc applying the sflats.


_________________________________________________________________________________________________

 

Part 2: World Coordinate Systems (WCS)
mscsetwsc

mscgetcatalog

msczero (and/or msccmatch)

msctvmark


All of the tasks to be run in this second part should be run on each object individulally and separately.  Meaning, do not do obj*.fits; specify and run each independently – run the task on obj1034.fits, run the task again on the next object obj1035.fits.

 


First you need to locate and get a naval stars database.  My advisor has one on his computer so I just copied it into the folder I was working in. 
                cp /d/glacier/ewilcots/GROUPS/WIYN/Vk1003_980909_36inch.db .

mscsetwcs:
Specify your database and specify the first object you will be working on in images.  As an example I will pretend I am working on processing an object called obj1011.fits.  this task will help set the wcs for the image to be used.
                PACKAGE = mscred
                   TASK = mscsetwcs

                images  =      obj1011.fits  Mosaic images
                database= Vk1003_980909_36inch.db  WCS database
                (ra     =                      telra) Right ascension keyword (hours)
                (dec    =               teldec) Declination keyword (degrees)
                (equinox=       telequin) Epoch keyword (years)
                (ra_offs=                     0.) RA offset (arcsec)
                (dec_off=                   0.) Dec offset (arcsec)
                (extlist=                         )
                (mode   =                   ql)


mscgetcatalog:

Next, run mscgetcatalog on the object to generate a coordinate catalog to be used in msczero.
                PACKAGE = mscred
                   TASK = mscgetcatalog

                input   =               obj1011.fits  List of Mosaic files
                output  =               cord1011  Output file of sources
                (magmin =                                5 ) Minimum magnitude
                (magmax =                             17) Maximum magnitude
                (catalog=      NOAO:USNO-A2) Catalog
                (rmin   =                                   30.) Minimum radius (arcmin)
                (mode   =                                   ql)


msczero:
We need to check that the image has the right world coordinate system.  We will run msczero on each amplifier of each object – this is unfortunately the longest and worst part of the process.  You’ll basically display one of the 8 parts of the image at a time, and then display the coordinate system for that piece.  A bunch of red circles will be displayed and you will place one of the red circles on the star it is suppose to be on.  You will then redisplay the coordinate circles and they should all now match up on the correct star.  At one point I ran into some sort of display problem, to fix it I typed…
                > set disable_wcs_maps=""
                > flpr
As I said, this is a one at a time process for best accuracy.  For input, specify your object and then what amplifier you are on (obj1011.fits[1], ... , obj1011.fits[8]).  For coords, specify the coordinate system you just generated.


                PACKAGE = mscred
                   TASK = msczero

                input   =            obj1011.fits[1]  List of mosaic exposures
                (extname=                     ) Extension name pattern
                (nframes=                    2) Number of frames to use
                (cbox   =                     11) Centering box size (see imcntr)
                (mark   =                  yes) Mark display?
                (logfile=            default) Log file for measurements
                                # MSCTVMARK Parameters
                coords  =               cord1011  List of coordinates
                (fields =                  1,2,3) Fields for RA, DEC, and ID
                (wcs    =                world) Coordinate type (logical|physical|world)
                (mtype  =               circle) Mark type
                radii   =                   20  Radii of concentric circles
                color   =                  204  Gray level of marks to be drawn
                label   =                   no  Label the marked coordinates
                (nxoffse=                   20) X offset in display pixels of number
                (nyoffse=                     0) Y offset in display pixels of number
                (pointsi=                      3) Size of mark type point in display pixels
                (txsize =                       2) Size of text and numbers in font units
                                # Task query and internal parameters
                ra      =                       RA (hours)
                dec     =                     DEC (degrees)
                id      =                       Identification
                mag     =                    Magnitude limit
                update  =                  yes  Update WCS zero point?
                updcoord=                yes  Update coordinate file?
                (fd1    =                          )
                (fd2    =      uparm$mscdisp1)
                (mode   =                   ql)

Run the task.

                >msczero

Amplifier will be displayed.  Once it is displayed click in the ds9 window and hit the letter m (for ‘mark’).  In terminal window you will be prompted to enter some information, but if you entered your parameters as I have, you’ll just have to hit [enter] four times.  The red circles should now be displayed.  Figure out what circles match up with which stars.  In the ds9 window select a red circle, hit ‘s’, select its star, hit z.  This will log the RA and DEC of the stars, so you should just have to hit [enter] twice.  The red circle should now be on/circling its star.  In the ds9 window hit ‘m’ again.  Hit [enter] again four times.  Now, all of the circles should be on their proper stars.  If they are not, you didn’t figure out what circle went with which star correctly.  When finished, hit ‘q’ in the ds9 window. 

Run msczero on each amplifier for each object.

 

When that is done, you can also run msccmatch.  Msccmatch is an interactive way to make sure the WCS is matched up correctly.  It can be run in addition to msczero or instead of msczero if your WSC somehow magically lined up without performing msczero.

                Note: I typically just ran msczero.  For most of my objects I didn’t find msccmatch to help much.

                PACKAGE = mscred
                   TASK = msccmatch

                input   =               obj1011.fits[1]  List of input mosaic exposures
                coords  =               cord1011  Coordinate file (ra/dec)
                (outcoor=             cordup1011) List of updated coordinate files
                (usebpm =                yes) Use bad pixel masks?
                (verbose=                 yes) Verbose?
                                                # Coarse Search
                (nsearch=                    60) Maximum number of positions to use in search
                (search =                    60.) Translation search radius (arcsec)
                (rsearch=                   0.2) Rotation search radius (deg)
                                # Fine Centroiding
                (cbox   =                      11) Centering box (pixels)
                (maxshif=                  3.5) Maximum centering shift to accept (arcsec)
                (csig   =                      0.1) Maximum centering uncertainty to accept (arcsec)
                (cfrac  =                     0.5) Minimum fraction of accepted centers
                (listcoo=                   yes) List centered coordinates in verbose mode?
                               # WCS Fitting
                (nfit   =                       50) Min for fit (>0) or max not found (<=0)
                (rms    =                    0.5) Maximum fit RMS to accept (arcsec)
                (fitgeom=              general) Fitting geometry
                (reject =                      3.) Fitting rejection limit (sigma)
                (update =                 yes) Update coordinate systems?
                (interac=                  yes) Interactive?
                (fit    =                       yes) Interactive fitting?
                (graphic=        stdgraph) Graphics device
                (cursor =                         ) Graphics cursor
                accept  =                  yes  Accept solution?
                (mode   =                   ql)
               
Run the task. 

                >msccmatch

The idea in this task is you will be removing points that do not fit well in the spread and finding a new order number that matches the points better.

An interactive window will pop up.  When this happens hit ‘x’.  Delete any points that don’t fit the pattern or are too far out with ‘d’.  Hit ‘f’ to apply.  Type ':xxorder *'  where * is a number typically 3, 4, or 5  whichever creates the best fit, remember to hit ‘f’ to apply.  Once you found a good fit move on; type ':yyorder *' where * is a number typically 3, 4, or 5  whichever creates the best fit, remember to hit ‘f’ to apply.  Hit ‘q’ to quit.  Then answer yes to save the changes you’ve made.  After each run, you’ll need to delete the cordup1011 file it creates.



msctvmark:
Once that is done, display your object in ds9.  We’re going to see how well you’re new WCS matches up with the image as a whole.  Running this task will display the red circles over your whole image.  Remember to specify your coordinate system for coords and make sure you display your image and the cords in the same ds9 frame.

                PACKAGE = mscred
                   TASK = msctvmark

                coords  =               cord1011  List of coordinates
                frame   =                    1  Display frame
                (output =                      ) Output file of pixel coordinates and labels
                (fields =                1,2,3) Fields for RA, DEC, and ID
                (wcs    =              world) Coordinate type (logical|physical|world)
                (mark   =              circle) Mark type
                (radii  =                     10) Radii of concentric circles
                (lengths=                    0) Lengths and width of concentric rectangles
                (font   =              raster) Default font
                (color  =                  204) Gray level of marks to be drawn
                (label  =                    no) Label the marked coordinates
                (nxoffse=                    0) X offset in display pixels of number
                (nyoffse=                    0) Y offset in display pixels of number
                (pointsi=                     3) Size of mark type point in display pixels
                (txsize =                      1) Size of text and numbers in font units
                (mode   =                   ql)

Does it look ok? If so, great! You’re done with this part! 

Stars and circles still not matching up?  Looks like you’re going to have to go back and try again…
_________________________________________________________________________________________________

Part 3: Getting the Combined Image
craverage

fixpix

                crmedian

                crgrow

                fixpix

mscimage

mscskysub

mscgetcatalog

mscimatch

mscstack



For part 3, I have already grouped my images to be combined into its own folder.  We’ll say I’m creating my r-band image of NGC383.  In my current folder I have the r-band images of NGC383 obj1011.fits, obj1012.fits, and obj1013.fits and of course a folder called ‘wscbackup’ containing  a copy of the images processed through part 1 and 2 – so that incase something goes wrong in combining I wont have to do part 2 again!

 

What we are going to do first is create some cosmic ray masks for our images.  I have not mastered the cosmic ray masks, but doing the following process got rid of enough of them.

 

 

craverage:

get to the package crutil (the location of the craverage command)

                >mscred

                >imred

                >crutil
First we need to create a script to use in craverage.  One way is to manually create one in your favorite text editor.  This is usually the way I did it – I like gedit so I went into a new word document and created a file called ‘craverage’ with the contents looking like this:

                craverage obj1011.fits[im1]  output="" crmask="crmask1011_1"
                craverage obj1011.fits[im2]  output="" crmask=" crmask1011_2"
                craverage obj1011.fits[im3]  output="" crmask=" crmask1011_3"
                craverage obj1011.fits[im4]  output="" crmask=" crmask1011_4"
                craverage obj1011.fits[im5]  output="" crmask=" crmask1011_5"
                craverage obj1011.fits[im6]  output="" crmask=" crmask1011_6"
                craverage obj1011.fits[im7]  output="" crmask=" crmask1011_7"
                craverage obj1011.fits[im8]  output="" crmask=" crmask1011_8"

                craverage obj1012.fits[im1]  output="" crmask="crmask1012_1"
                craverage obj1012.fits[im2]  output="" crmask=" crmask1012_2"
                craverage obj1012.fits[im3]  output="" crmask=" crmask1012_3"
                craverage obj1012.fits[im4]  output="" crmask=" crmask1012_4"
                craverage obj1012.fits[im5]  output="" crmask=" crmask1012_5"
                craverage obj1012.fits[im6]  output="" crmask=" crmask1012_6"
                craverage obj1012.fits[im7]  output="" crmask=" crmask1012_7"
                craverage obj1012.fits[im8]  output="" crmask=" crmask1012_8"

                craverage obj1013.fits[im1]  output="" crmask="crmask1013_1"
                craverage obj1013.fits[im2]  output="" crmask=" crmask1013_2"
                craverage obj1013.fits[im3]  output="" crmask=" crmask1013_3"
                craverage obj1013.fits[im4]  output="" crmask=" crmask1013_4"
                craverage obj1013.fits[im5]  output="" crmask=" crmask1013_5"
                craverage obj1013.fits[im6]  output="" crmask=" crmask1013_6"
                craverage obj1013.fits[im7]  output="" crmask=" crmask1013_7"
                craverage obj1013.fits[im8]  output="" crmask=" crmask1013_8"

Make sure that there is a command line for each object for each of the 8 parts you have in your folder.

 

Another way to get this is to use the task mscstat.  First, get mscstat parameters look like this (but do not run yet):
                PACKAGE = mscred
                   TASK = mscstat

                images  =            obj*.fits  Images
                (extname=                       ) Extension name selection
                (gmode  =                   no) Global mode statistics?
                (fields =                image) Fields to be printed
                (lower  =                INDEF) Lower cutoff for pixel values
                (upper  =                INDEF) Upper cutoff for pixel values
                (binwidt=                  0.1) Bin width of histogram in sigma
                (format =                  yes) Format output and print column labels?
                (fd1    =                            )
                (fd2    =                            )
                (mode   =                     ql)
enter the command…      

                mscstat > cravg_script
This will make the start of script, then double check in a editor that it looks like the one shown above.

 

Now that you have your script ready you can run craverage.
                PACKAGE = crutil
                   TASK = craverage

                input   =           List of input images
                output  =           List of output images
                (crmask =          ) List of output cosmic ray and object masks
                (average=         ) List of output block average filtered images
                (sigma  =           ) List of output sigma images
                (navg   =         7) Block average box size
                (nrej   =          3) Number of high pixels to reject from the average
                (nbkg   =        5) Background annulus width
                (nsig   =        50) Box size for sigma calculation
                (var0   =        0.) Variance coefficient for DN^0 term
                (var1   =        0.) Variance coefficient for DN^1 term
                (var2   =        0.) Variance coefficient for DN^2 term
                (crval  =         1) Mask value for cosmic rays
                (lcrsig =      10.) Low cosmic ray sigma outside object
                (hcrsig =  3.75) High cosmic ray sigma outside object
                (crgrow =     0.) Cosmic ray grow radius     
                (objval =        0) Mask value for objects
                (lobjsig=     10.) Low object detection sigma
                (hobjsig=    5.5) High object detection sigma
                (objgrow=    0.) Object grow radius
                (mode   =      ql)

enter the command    

                cl < cravg_script
You are now creating cosmic ray masks.  In my notes I originally wrote “this takes forever and a year...”  so be prepared to let your computer run for a bit.  Sometimes IRAF kicked out a floating point error, would stop creating the masks… but I simply ran the task again and it would work just fine starting right where it left off.

 

fixpix:
Next we need two more text files.  One called fixpix_in that has all the objects and one called crmask that lists all the crmasks that have been created.  Again, you can create the text files by hand in a text editor or using mscstat.  If using mscstat, enter the command…  
                mscstat > fixpix_in
For this example folder it looks like this:
                obj1011.fits[im1]
                obj1011.fits[im2]
                obj1011.fits[im3]
                obj1011.fits[im4]
                obj1011.fits[im5]
                obj1011.fits[im6]
                obj1011.fits[im7]
                obj1011.fits[im8]
                obj1012.fits[im1]
                obj1012.fits[im2]
                obj1012.fits[im3]
                obj1012.fits[im4]
                obj1012.fits[im5]
                obj1012.fits[im6]
                obj1012.fits[im7]
                obj1012.fits[im8]
                obj1013.fits[im1]
                obj1013.fits[im2]
                obj1013.fits[im3]
                obj1013.fits[im4]
                obj1013.fits[im5]
                obj1013.fits[im6]
                obj1013.fits[im7]
                obj1013.fits[im8]
Again, you can create the text file for crmask by hand in a text editor or using mscstat.  If using mscstat, enter the command…       
                mscstat > crmask

The text file for crmask should look like this…

                crmask1011_1.fits[im1]
                crmask1011_2.fits[im2]
                crmask1011_3.fits[im3]
                crmask1011_4.fits[im4]
                crmask1011_5.fits[im5]
                crmask1011_6.fits[im6]
                crmask1011_7.fits[im7]
                crmask1011_8.fits[im8]
                crmask1012_1.fits[im1]
                crmask1012_2.fits[im2]
                crmask1012_3.fits[im3]
                crmask1012_4.fits[im4]
                crmask1012_5.fits[im5]
                crmask1012_6.fits[im6]
                crmask1012_7.fits[im7]
                crmask1012_8.fits[im8]
                crmask1013_1.fits[im1]
                crmask1013_2.fits[im2]
                crmask1013_3.fits[im3]
                crmask1013_4.fits[im4]
                crmask1013_5.fits[im5]
                crmask1013_6.fits[im6]
                crmask1013_7.fits[im7]
                crmask1013_8.fits[im8]
make sure the objects and crmasks orderings match up. Meaning, obj1011[1] is first and so is crmask1011[1]
Run the task fixpix with the following parameters.  This will apply the crmasks to your image.
                PACKAGE = proto
                   TASK = fixpix

                images  =           @fixpix_in  List of images to be fixed
                masks   =              @crmask  List of bad pixel masks
                (linterp=                 INDEF) Mask values for line interpolation
                (cinterp=                INDEF) Mask values for column interpolation
                (verbose=                   yes) Verbose output?
                (pixels =                       no) List pixels?
                (mode   =                      ql)

If you held onto to the original images, you can now compare the before the crmasks were applied to the image with the crmasks applied.

 

crmedian:

There is another way or an additional way to get rid of some cosmic rays.  The other is using crmedian instead of craverage as well as applying a growing of what is to be covered.

 

made a script like that for craverage except for crmedian
                crmedian obj1011.fits[im1]  output="" crmask="crmask1011_1"
                crmedian obj1011.fits[im2]  output="" crmask=" crmask1011_2"
                crmedian obj1011.fits[im3]  output="" crmask=" crmask1011_3"
                crmedian obj1011.fits[im4]  output="" crmask=" crmask1011_4"
                crmedian obj1011.fits[im5]  output="" crmask=" crmask1011_5"
                crmedian obj1011.fits[im6]  output="" crmask=" crmask1011_6"
                crmedian obj1011.fits[im7]  output="" crmask=" crmask1011_7"
                crmedian obj1011.fits[im8]  output="" crmask=" crmask1011_8"

                crmedian obj1012.fits[im1]  output="" crmask="crmask1012_1"
                crmedian obj1012.fits[im2]  output="" crmask=" crmask1012_2"
                crmedian obj1012.fits[im3]  output="" crmask=" crmask1012_3"
                crmedian obj1012.fits[im4]  output="" crmask=" crmask1012_4"
                crmedian obj1012.fits[im5]  output="" crmask=" crmask1012_5"
                crmedian obj1012.fits[im6]  output="" crmask=" crmask1012_6"
                crmedian obj1012.fits[im7]  output="" crmask=" crmask1012_7"
                crmedian obj1012.fits[im8]  output="" crmask=" crmask1012_8"

                crmedian obj1013.fits[im1]  output="" crmask="crmask1013_1"
                crmedian obj1013.fits[im2]  output="" crmask=" crmask1013_2"
                crmedian obj1013.fits[im3]  output="" crmask=" crmask1013_3"
                crmedian obj1013.fits[im4]  output="" crmask=" crmask1013_4"
                crmedian obj1013.fits[im5]  output="" crmask=" crmask1013_5"
                crmedian obj1013.fits[im6]  output="" crmask=" crmask1013_6"
                crmedian obj1013.fits[im7]  output="" crmask=" crmask1013_7"
                crmedian obj1013.fits[im8]  output="" crmask=" crmask1013_8"

Make sure that there is a command line for each object for each of the 8 parts you have in your folder.

                PACKAGE = crutil
                   TASK = crmedian

                input   =                       Input image
                output  =                       Output image
                (crmask =                      ) Output cosmic ray mask
                (median =                     ) Output median image
                (sigma  =                       ) Output sigma image
                (residua=                      ) Output residual image
                (var0   =                     0.) Variance coefficient for DN^0 term
                (var1   =                     0.) Variance coefficient for DN^1 term
                (var2   =                     0.) Variance coefficient for DN^2 term
                (lsigma =                 10.) Low clipping sigma factor
                (hsigma =                4.5) High clipping sigma factor
                (ncmed  =                   5) Column box size for median level calculation
                (nlmed  =                    5) Line box size for median level calculation
                (ncsig  =                    25) Column box size for sigma calculation
                (nlsig  =                     25) Line box size for sigma calculation
                (mode   =                   ql)
then run crmedian with command
                cl < crmed_script

 

 

crgrow:
Next, run crgrow on all of the crmasks that were made
                PACKAGE = crutil
                   TASK = crgrow

                input   =       crmask1011_8[1]  Input cosmic ray masks
                output  =       crmask1011_8[1]  Output cosmic ray masks
                (radius =                       1.) Grow radius
                (inval  =                INDEF) Input mask value to grow (INDEF for any)
                (outval =              INDEF) Output grown mask value (INDEF for any)
                (mode   =                    ql)

 

fixpix:
Apply with fixpix task.  Remember to make sure fixpix_in and crmask scripts match up correctly.
                PACKAGE = proto
                   TASK = fixpix

                images  =           @fixpix_in  List of images to be fixed
                masks   =              @crmask  List of bad pixel masks
                (linterp=                 INDEF) Mask values for line interpolation
                (cinterp=                INDEF) Mask values for column interpolation
                (verbose=                   yes) Verbose output?
                (pixels =                       no) List pixels?
                (mode   =                      ql)

If you held onto to the original images, you can now compare the before the crmasks were applied to the image with the crmasks applied.  Also might be a good idea to check which type of cosmic ray masks works best for your data.  I applied both craverage and crmedian with crgrow to my images.


mscimage:

run mscimage on each object.  This step makes it so IRAF doesn’t identify the image as 8 separate amplifiers anymore.

                PACKAGE = mscred
                   TASK = mscimage
               
                input   =               obj1011  List of input mosaic exposures
                output  =               mos1011  List of output images
                (format =                image) Output format (image|mef)
                (pixmask=                   yes) Create pixel mask?
                (verbose=                   yes) Verbose output?
                                # Output WCS parameters
                (wcssour=              image) Output WCS source (image|parameters)
                (referen=                          ) Reference image
                (ra     =                   INDEF) RA of tangent point (hours)
                (dec    =                 INDEF) DEC of tangent point (degrees)
                (scale  =                 INDEF) Scale (arcsec/pixel)
                (rotatio=               INDEF) Rotation of DEC from N to E (degrees)
                              # Resampling parmeters
                (blank  =                       0.) Blank value
                (interpo=              sinc17) Interpolant for data
                (minterp=              linear) Interpolant for mask
                (boundar=       constant) Boundary extension
                (constan=                     0.) Constant boundary extension value
                (fluxcon=                    no) Preserve flux per unit area?
                (ntrim  =                        7) Edge trim in each extension
                (nxblock=                4200) X dimension of working block size in pixels
                (nyblock=                2100) Y dimension of working block size in pixels
                                 # Geometric mapping parameters
                (interac=                    no) Fit mapping interactively?
                (nx     =                        10) Number of x grid points
                (ny     =                        20) Number of y grid points
                (fitgeom=          general) Fitting geometry
                (xxorder=                      4) Order of x fit in x
                (xyorder=                      4) Order of x fit in y
                (xxterms=                 half) X fit cross terms type
                (yxorder=                      4) Order of y fit in x
                (yyorder=                      4) Order of y fit in y
                (yxterms=                 half) Y fit cross terms type
                (fd_in  =                           )
                (fd_ext =                          )
                (fd_coor=                        )
                (mode   =                     ql)
               

mscskysub:

Make a text file called inskysub.  File should list all of the mos*.fits files.

                mos1011.fits

                mos1012.fits

                mos1013.fits
make a text file called outskysub. List new names for the files like follows…
                mos1011_2.fits
                mos1012_2.fits
                mos1013_2.fits

Once you’ve made the files, run mscskysub

                PACKAGE = mscred
                   TASK = mscskysub

                input   =            @inskysub  Input images to be fit
                output  =           @outskysub  Output images
                xorder  =                    2  Order of function in x
                yorder  =                    2  Order of function in y
                (type_ou=             residual) Type of output (fit,residual,response,clean)
                (functio=                        leg) Function to be fit (legendre,chebyshev,spline3)
                (cross_t=                       yes) Include cross-terms for polynomials?
                (xmedian=                    100) X length of median box
                (ymedian=                    100) Y length of median box
                (median_=                     50.) Minimum fraction of pixels in median box
                (lower  =                          0.) Lower limit for residuals
                (upper  =                          0.) Upper limit for residuals
                (ngrow  =                          0) Radius of region growing circle
                (niter  =                             0) Maximum number of rejection cycles
                (regions=                   mask) Good regions (all,rows,columns,border,sections,c
                (rows   =                            *) Rows to be fit
                (columns=                        *) Columns to be fit
                (border =                        50) Width of border to be fit
                (section=                            ) File name for sections list
                (circle =                               ) Circle specifications
                (mask   =                  BPM) Mask
                (div_min=               INDEF) Division minimum for response output
                (mode   =                        ql)


mscgetcatalog:
run mscgetcatalog with following parameters on only one/the first mos*.fits image so that the one image can be used as a reference for the rest.
                PACKAGE = mscred
                   TASK = mscgetcatalog

                input   =         mos1011_2.fits  List of Mosaic files
                output  =       mos1011_2.usnoA2  Output file of sources
                (magmin =                   17.) Minimum magnitude
                (magmax =                  22.) Maximum magnitude
                (catalog=         NOAO:USNO-A2) Catalog
                (rmin   =                       21.) Minimum radius (arcmin)
                (mode   =                       ql)


mscimatch: 

This is run on all the files, it will apply coordinates to all the mos*_2.fits files so all images are using the same reference.
                PACKAGE = mscred
                   TASK = mscimatch

                input   =           @outskysub  List of images
                coords  =      mos1011_2.usnoA2  File of coordinates
                (bpm    =                BPM) List of bad pixel masks
                (measure=                     ) Measurment file
                (scale  =                   yes) Determine scale?
                (zero   =                     no) Determine zero offset?
                (box1   =                    21) Inner box size for statistics
                (box2   =                    51) Outer box size for statistics
                (lower  =               -100.) Lower limit for good data
                (upper  =             INDEF) Upper limit for good data
                (niterat=                      3) Number of sigma clipping iterations
                (sigma  =                    3.) Sigma clipping factor
                (interac=                   no) Interactive?
                (verbose=                yes) Verbose?
                accept  =                  yes  Accept scaling and update images?
                (mode   =                   ql)


mscstack:

Hurray! The last step! This step will stack all the images together to give you your final combined image!   Name the output something that will be easy to identify; I chose what is being observed and what filter its in.
                PACKAGE = mscred
                   TASK = mscstack

                input   =           @outskysub  List of images to combine
                output  =              NGC383_r  Output image
                (headers=                      ) List of header files (optional)
                (bpmasks=           *_bpm.fits) List of bad pixel masks (optional)
                (rejmask=                      ) List of rejection masks (optional)
                (nrejmas=                      ) List of number rejected masks (optional)
                (expmask=                     ) List of exposure masks (optional)
                (sigmas =                        ) List of sigma images (optional)
                (combine=        median) Type of combine operation (median|average)
                (reject =              ccdclip) Type of rejection
                (masktyp=   goodvalue) Mask type
                (maskval=                    0.) Mask value
                (blank  =                       0.) Value if there are no pixels
                (scale  =         !mscscale) Image scaling
                (zero   =          !msczero) Image zero point offset
                (weight =               none) Image weights
                (statsec=                         ) Image section for computing statistics
                (lthresh=                      1.) Lower threshold
                (hthresh=             INDEF) Upper threshold
                (nlow   =                        1) minmax: Number of low pixels to reject
                (nhigh  =                        1) minmax: Number of high pixels to reject
                (nkeep  =                       1) Minimum to keep (pos) or maximum to reject (neg)
                (mclip  =                    yes) Use median in sigma clipping algorithms?
                (lsigma =                      5.) Lower sigma clipping factor
                (hsigma =                   2.9) Upper sigma clipping factor
                (rdnoise=           rdnoise) ccdclip: CCD readout noise (electrons)
                (gain   =                     gain) ccdclip: CCD gain (electrons/DN)
                (snoise =                       0.) ccdclip: Sensitivity noise (fraction)
                (sigscal=                     0.1) Tolerance for sigma clipping scaling corrections
                (pclip  =                    -0.5) pclip: Percentile clipping parameter
                (grow   =                      0.) Radius (pixels) for neighbor rejection
                (mode   =                     ql)


You’re done! Check out the image!



Comments