Adding a WCS to your fits image

So, adding a WCS is not actually all that difficult. But in my experience finding the process clearly explained is more challenging. I finally found a process that works well for me (with thanks going to Kyle Johnston, U. Waterloo, for outlining the process). Also, there are probably multiple ways to do this, but here is my attempt to lay it out clearly. I hope it helps!

(for a more formatted version, try this pdf file)

Steps for Adding a WCS to your FITS image:

1) Load your image into ds9

2) Create a new frame and load in a DSS image of your field (Analysis --> Image Servers --> SAO DSS (or whichever server you prefer))

It helps to tile the frames and arrange their zooms, etc. so that they look as identical as possible.

3) On the DSS frame, load in reference stars (usually USNO; Analysis --> Catalogs --> USNO-A2.0 Catalog).

It will automatically plot small circles around every USNO object in the field.

4) Change the pointer to "catalog" (Edit --> Catalog). This enables you to click to select each circle.

5) start "imexam" at the iraf prompt

6) Open up a new text file

7) Search for stars that overlap between the USNO catalog and your frame. Also, make sure they are observable (i.e. not saturated or too faint) on both frames.

8) Select one such star on your frame; mouse over and type , to get iraf to print the pixel coordinates of the center of that star to the screen. Copy those coordinates into the text file.

9) Click on the corresponding star in the DSS frame, making sure the circle of the USNO star is selected. This will highlight the appropriate row in the USNO table that popped up when you loaded the stars in ds9. Bring the table to the front and copy the text in that row. Paste it into the text file, on the same line as the x, y coordinates from imexam. You may need to delete extraneous text to just leave the RA & dec.

10) Repeat for as many stars as possible. It's good to get 20 or even more, if possible. Also, try to get the best coverage across the frame, as this will improve the quality of your solution.

At this point, you should have a nice text file with x, y, RA, Dec in it. This is your input into ccmap, the iraf task to find the WCS solution. Note that you might need to convert RA from degrees to hours (just multiply by 24/360, or 0.08)

the iraf task you need, ccmap is found in images.imcoords

11) In iraf, edit the ccmap parameters as follows (obviously, change the file names to reflect your own)

input = "starlist.dat" The input coordinate files
database = "n3311_coordfit.db" The output database file
(solutions = "") The database plate solution names
(images = "imask") The input images
(results = "results_summary") The optional results summary files
(xcolumn = 1) Column containing the x coordinate
(ycolumn = 2) Column containing the y coordinate
(lngcolumn = 3) Column containing the ra / longitude
(latcolumn = 4) Column containing the dec / latitude
(xmin = INDEF) Minimum logical x pixel value
(xmax = INDEF) Maximum logical x pixel value
(ymin = INDEF) Minimum logical y pixel value
(ymax = INDEF) Maximum logical y pixel value
(lngunits = "hours") Input ra / longitude units
(latunits = "degrees") Input dec / latitude units
(insystem = "j2000") Input celestial coordinate system
(refpoint = "coords") Source of the reference point definition
(lngref = "INDEF") Reference point ra / longitude telescope coordinate
(latref = "INDEF") Reference point dec / latitude telescope coordinate
(refsystem = "INDEF") Reference point telescope coordinate system
(lngrefunits = "INDEF") Reference point ra / longitude units
(latrefunits = "INDEF") Reference point dec / latitude units
(projection = "tan") Sky projection geometry
(fitgeometry = "general") Fitting geometry
(function = "polynomial") Surface type
(xxorder = 2) Order of xi fit in x
(xyorder = 2) Order of xi fit in y
(xxterms = "half") Xi fit cross terms type
(yxorder = 2) Order of eta fit in x
(yyorder = 2) Order of eta fit in y
(yxterms = "half") Eta fit cross terms type
(maxiter = 0) The maximum number of rejection iterations
(reject = 3.) Rejection limit in sigma units
(update = no) Update the image world coordinate system ?
(pixsystem = "logical") Input pixel coordinate system
(verbose = yes) Print messages about progress of task ?
(interactive = yes) Fit the transformation interactively ?
(graphics = "stdgraph") Default graphics device
(cursor = "") Graphics cursor
(mode = "ql")

This will pop up a window with the data points. You can interactively examine your fit by typing "x" or "y" to see residuals of the fit. These should be small (a few tenths of an arcsec at most, depending on the quality you are aiming for). You can mouse over bad data points and click "d" to delete them. Then type "f" to refit the data without those points. You can see a complete list of command options by typing "?"

NOTE: ccmap does NOT actually apply the WCS solution, it just calculates it!

12) Apply the newly found WCS to your image using ccsetwcs in iraf. The parameters look something like this:

images = "imask" The input images
database = "n3311_coordfit.db" The input database file
solutions = "imask" The input plate solutions
(xref = INDEF) The x reference pixel
(yref = INDEF) The y reference pixel
(xmag = INDEF) The x axis scale in arcsec / pixel
(ymag = INDEF) The y axis scale in arcsec / pixel
(xrotation = INDEF) The x axis rotation angle in degrees
(yrotation = INDEF) The y axis rotation angle in degrees
(lngref = INDEF) The ra / longitude reference coordinate in lngunits
(latref = INDEF) The dec / latitude reference coordinate in latunits
(lngunits = "") The ra / longitude reference coordinate units
(latunits = "") The dec / latitude reference coordinate units
(transpose = no) Transpose the computed image wcs ?
(projection = "tan") The sky projection geometry
(coosystem = "j2000") The celestial coordinate system
(update = yes) Update the image world coordinate system ?
(pixsystem = "logical") The input pixel coordinate system
(verbose = yes) Print messages about actions taken by the task ?
(mode = "ql")

Voila! You should have a nice WCS now attached to your image.

13) Check the WCS by loading your image (with the new WCS) into ds9. Then, load in the USNO stars as in step 3, except on your own image this time, instead of the DSS image. If all went well, the circles should appear nicely over their corresponding stars.

Permalink • Posted in: astronomyComments (1)


scige Jan 29, 2009

Thanks for the Human-Readable way you wrote this... =)