automatch:
Finding Precise Coordinates for Images

John Thorstensen
Last revised: July 2001

Introduction

This document describes software installed on agung and hill (the observer's Linux workstations at MDM) which sets up the so-called world coordinate system (or WCS) in your direct image headers. The WCS is a precise RA and dec solution in the image header. Once an image has the WCS set up correctly, you can instantly get a precise RA and dec for any point on your image. A rich set of applications use the WCS for mosaicking, aligning optical data with images from other sources, and so on.

Briefly, the procedure described here automatically (1) finds stars in your pictures (2) matches them up to entries in the huge USNO A2.0 catalog of stars, and (3) finally puts the coordinate solutions from these stars into the headers. The software is highly automated and can process many images quickly and conveniently, or can be used on single images. The automatic algorithms don't always work perfectly, but the software includes tools which can ferret out these cases and help you fix them efficiently.

The procedures here come in two flavors - single-image or batch. Doing many images in batch is efficient and also makes possible some useful consistency checks.

The first part of this document gives important background information and quick instructions. The rest gives more details, including suggestions for how to proceed when things don't work (this material is mostly toward the end).

Important Considerations - READ THIS

• These procedures don't yet work on 8k images - though you may have some success with individual images from the 8k with some hand massaging. (8k images have no telescope coordinates in header, and the software doesn't have a map of where the individual chips are.)

• In order to find stars automatically, the telescope coordinates in your image headers have to be reasonably accurate (like 1 arcmin).

• Your images must be in FITS format with .fits filename extension. Furthermore, there cannot be any identically-named images with the .imh filename extensions (those are obsolete anyway!), or IRAF gets very confused.

• Images should have the final geometry (cropping), but they need not be flatfielded.

• The procedures run from within IRAF, though most of the work is done outside IRAF using shell escapes (!).

• The detector pixels must be square and not rectangular (true of all modern CCDs, I think).

• If you're processing a bunch of images from several different setups, group them in similar bunches (same chip, telescope, and exact rotation) for consistency checks.

• The image scale must be known reasonably accurately (the header values are generally good enough), but there's no restriction as to parity or rotation. If the image scale isn't known, you'll have to do the first image by hand and use hedit to insert the image scale in the headers (see end of the document).

• Conventions used in this manual.
  This typeface indicates computer input and output.

  A * is sometimes used to substitute for something obvious, like an image name.

  I use some generic names like imagelist - this is a stand-in for a file of image names, which you can call whatever you want.

  Similarly, imagename is the (root) name of a generic image. However, when programs write scripts, the scripts always have the same name, so if you're running a machine-generated script, you just type literally from the instructions in this document. (Generally the programs end with a little note about what to do next, which should be helpful.)

Doing a Single Image

1. Copy the image you want to the linux box. I'll refer to it as imagename. As mentioned earlier the filename itself must have the .fits extension, but the procedures which follow should work whether or not you include the .fits in the command.

2. Get into IRAF. From within IRAF, type

       !make1match.p imagename
   

   followed by

       cl < matchscript1
   

   This will automatically do the following:
   1. Run a program to check the image header and, if necessary, preserve the original telescope coordinates, which are overwritten later;

   2. Generate a catalog of pixel coordinates of stars found in the image (using SExtractor and a program to filter its output);

   3. Look up stars in the vicinity in a locally-resident, abridged version of the USNO A2.0 catalog;

   4. attempt to establish a match between the stars detected in the picture and USNO catalog stars;

   5. If an apparent match is found, print out some information and ask you whether you want to accept the match. Note carefully how many stars were fit (compared to how many were there), the rotation angle, parity, and scale. These last should be stable during a run unless you play with the rotator.

   6. If you accept the match, output files are written with the appropriate information, and a script is produced which will load the information in the header (next step).

3. If the match solution looks to be good, load the WCS in the header by typing

      cl < wcs_script1
   

   Having this as a separate step offers an extra measure of safety against spurious matches. If you imhead imagename long+, you should now see various WCS keywords at the end of the image header, stuff like CTYPE1, CRVAL1, CRPIX1, CD1_1, and so on. If not there was some kind of trouble writing to the image.

4. You're now done! But, you may wish to check the match. The match is almost certainly good if (a) the program claims that there were a fair number of USNO stars in the field (like a dozen or so) and most (e.g., 8) were found and fitted; (b) the rms residual is small (like 0.5 arcsec); (c) the rotation angle is as expected (generally close to a multiple of 90 degrees, and nearly equal to other images taken with the same setup); (d) the image scale is equal to a few parts per thousand with other images taken with the same setup. Major violations of these conditions usually mean the match is wrong, or that the configuration of stars is particularly unfavorable.

   Sometimes you won't be certain, and if you're doing something critical (e.g. reporting coordinates of your hot new object for the IAUC) you'll want to check it anyway. Here's how to check an existing solution:

   • Display the image from IRAF in an xterm or DS9 window - display imagename 1

   • In an xterminal (not in IRAF), type

         handmatchsno imagename
     

     This will ask some questions, then pop a chart of the USNO stars in the area of your picture.

   • Type g to get the solution information (from the imagename.radec.match.wcs file, and run the solution.

   • Type o to overplot little boxes and XY coordinates on the USNO plot.

   • Back in IRAF, type tvmark 1 xytemp, which will mark the same stars on the plot. From this (and the fit information) you should get a pretty good idea of how reliable the solution is. handmatchusno has some options, described more fully at the end of this document, for entering new matches and editing out bad ones. If you're happy with the solution you have, type x to exit without overwriting anything; if you do want to save changes, type q to write out results and quit.

     Note that if you display the wcs'd image in DS9 (by reading it in directly, not in IRAF), the RA and dec of the cursor will display real-time! (If the frame appears all black, select scale, then zscale to reproduce the IRAF autoscaling. DS9 knows the whole range of pixel values, not just the 200 or so which ximtool uses.).

What it if fails?

The last few pages of this document for cleaning up cases in which the automatic matching failed. Often the header information is screwed up. Sometimes some hand work with SExtractor output will give a better catalog of detected stars to work with than the automated version does. Options available in handmatchusno can be used to rescue tough cases with a little effort on your part - as little as possible, I hope!

Doing a batch of images.

1. First, copy the images you want to do to the data disk on the Linux box.

2. Make a list of the images; I'll refer to it as imagelist but any name will do. There are several ways to do this, like files *.fits > imagelist or hsel *.fits $I 'naxis1 > 500' > imagelist. or you could use a text editor. Be careful not to leave a blank line at the end of the list - it'll generate some relatively harmless weirdness if you do. The program lop off any .fits extensions, so don't sweat those.

3. Type

       !makematch.p imagelist
   

   followed by

       cl < matchscript
   

   This does a great deal of work and may take tens of minutes to run if you have many images. The steps for each image are essentially as described earlier for single images. If even a plausible match is found (>3 stars), an output file is written without querying the user.

4. Were all your images taken at the same rotator angle? If so, type

       !checkwcs imagelist 1
   

   If the rotator angle was not consistent, omit the trailing 1. Then type

       cl < wcs_script > wcsout
   

   to load the wcs information in the headers.

5. The checkwcs step applies some consistency checks to the procedure. If you want, you can check the process further:

   • headaudit.tmp keeps a record of the step in which the headers are automatically queried and, if appropriate, the original telescope coordinates are saved.

   • unmatched_pix will be a list of pictures for which the procedure failed.

   • match_diagn.tmp gives details of the star matching results.

   • checkwcs also writes a file called wcsaudit which contains details of the solutions to be used in the WCS.

   • wcsout, if captured as suggested above, will have gory details of the WCS loading procedure.

6. The last few pages of this document have some hints for cleaning up cases in which the automatic matching failed. Often the header information is wrong or incomplete. Sometimes the algorithms which detect stars in images get confused, in which case a little hand work with SExtractor output will give a better catalog of detected stars to work with than the automated version does. Occasionally the automatic matching algorithm gets fooled into producing a false correspondence between detected stars and USNO catalog entries. A semi-automated program handmatchusno can be used to rescue tough cases. Just be sure not to re-run the keyword saving step (3) above!

7. When you're completely done, you can erase the files ending with .radec.match.wcs, .hsel, and .cat. The .db files are used by IRAF for some applications of the WCS - don't ask me why they don't just use the information in the images. So you should may want to keep 'em. [Note that the procedure affects some important header information. At MDM, the keywords RA, DEC, and EQUINOX are used by CCDCOM to record the telescope coordinates and their equinox (i.e., the numbers derived from the encoder readings). The WCS procedure commandeers the EQUINOX keyword for its own purposes, and there's no way around it except to reverse-engineer the WCS, which is a sizable job. For this reason the script runs a program which probes the headers and, if it appears appropriate, preserves `safe' copies of the original information in the new keywords TELRA, TELDEC, and TELEQNOX. This alone would leave the problem that the RA and DEC could disagree with the (new) EQUINOX keyword - if you'd taken your data with the telescope readout in 1950, say, the RA and DEC would be for 1950, but the EQUINOX would read 2000 (which is the WCS default). To avoid this problem, I've set up the wcs_script to load the computed coordinates of the exact field center into the RA and DEC fields.]

Gory Details

This section describes how these procedures work in more detail, which may be of interest to some. If the canned procedure doesn't work for you, you may be able to figure out the problem and fix it using this more detailed material.

Essentially, the procedures uses scripts or programs to write other scripts, which check all sorts of things and plug other software together like tinker toys. The über-scripts are partly generated using Perl (the ones with the .p suffix) or compiled C programs. These programs generate simple IRAF cl scripts, which you execute from within IRAF by typing `cl < scriptname'. If you look at the cl scripts, you'll find that they contain a fair number of IRAF shell escapes ! - a command preceded by this is simply issued to the usual user environment outside IRAF. I'll go through everything here in the order it happens.

Preparing the headers.

The aformentioned problem with the header keywords means that we have to save the telescope coordinate information separately in order to avoid destroying it. To do this I've written a C program prephead, which probes the headers and which generates an IRAF cl script to save the information. It only does this if it finds that the original information has not already been preserved (to avoid overwriting it) and that the WCS has not been set (in which case the pointing information is already lost if not yet preserved). The original telescope coordinates (born as RA, DEC, and EQUINOX) is copied into new keywords TELRA, TELDEC, and TELEQNOX in the following lines:

   hedit imagename teleqnox "(equinox)" add+ update+ delete-
   hedit imagename telra "(ra)" add+ update+ delete-
   hedit imagename teldec "(dec)" add+ update+ delete-

The header preparation script should only be run once. Therefore, when it is run, its last act is to erase itself!

What the Matchsript Does

The makematch.p perl script reads in the image list, strips off the .fits extensions if they're present, and writes out a script called matchscript, which executes the same commands for each image. When executed, the matchscript cl does the really heavy work, which I'll describe here step by step. I'll use imagename for the root image name (that is, minus the .fits) of the image being processed. [The single-image version make1match.p is very similar, but writes single-image versions of some of the commands, flagged by an extra argument s on the command line. This is all hidden from the user.]

The hsel command reads the following header information:

 
    telra  teldec  teleqnox  naxis1  naxis2  secpix1  filter

The sex command

Ah, those French, always witty ones. This runs the incredible SExtractor program, which chews through the picture and generates a list of objects found within it. There are a number of configuration files which I've tuned up for the task at hand and put in a centralized place (/usr/local/pkg/thor/sexdef). The configuration files are set up to: (1) Turn off the option which convolves the image with a filter before detecting objects, because we want to preserve the value of the central pixel. (2) sets the threshold for detection somewhat above the limit, to avoid spurious detections. (3) Sets the parameters which are to be reported for each detection in the output; a specific set is expected by later software. SExtractor generates a file called test.cat, which is transient. For more detail, see:

http://www.eso.org/science/eis/eis_doc/sex2/sex2html/sex2_doc.htm

The cat2seeing command

This reads the test.cat file generated by SExtractor and filters it to find real stars. It generates as output a file called imagename.cat. First, it filters the output for cosmic rays using three different tests (the `magic numbers' given here are defined constants in the code, which can be changed easily on recompilation. If you need to do this, copy and make your own private versions). (1) Objects with FWHM < 0.9 are rejected. (2) If the central pixel contains more than 40 per cent of the total flux, the object is rejected. This could cause the procedure to fail if the seeing gets so good that stars are undersampled. (3) If the image is badly out of round, the image is rejected. This will cause failure if the image is badly trailed. Another filter looks for grossly too many detections in any one single column, which is then rejected - this guards against spurious detections along bad columns. The FWHM seeing is then estimated from the mode of the FWHMa of the remaining detections, and only detections with FWHM within a reasonable margin of the modal FWHM are retained. This step can fail if there are very few bona-fide stars on an image. There's an exception: objects which are larger than the FWHM, and which have a maximum pixel with more than 30,000 counts, are taken to be saturated stars, and are accepted. If you've done something screwy with the middles of the bright stars (e.g., lopped off the peaks at some smaller data value), they won't be counted properly. Finally, detections near the limit are rejected. The output imagename.cat file contains the modal FWHM on the first line (which starts with a #), followed by the stars in order from the brightest to the faintest.

There are other procedures for finding stars in images (e.g., DAOFIND), but the ones I've seen need an initial guess as to the FWHM. This is fine when you're poring over a precious summary image of some cluster, but it's far too much work with hundreds of images. The SExtractor procedure used here usually does a reasonable job with no human intervention. DAOFIND is also problematic when it comes to saturated images, which the modified SExtractor command seems to deal with reasonably well.

The matchusno command

This takes as input the imagename, then does the following. (1) It grabs the imagename.hsel file generated earlier, and from that information infers the location, scale, and size of the image. (2) It goes to a compressed version of the USNO A2.0 star catalog (more on that later) and grabs the stars within a box which, to allow for pointing error, is 1 arcmin larger than the largest dimension of the image (so 30" larger at each edge). Using the filter (UBVRI) info from the imagename.hsel file, and the USNO b and r magnitudes, it computes rough magnitudes in the image's native passband, and then sorts the USNO stars in brightness order. (3) It reads in the imagename.cat file containing the filtered SExtractor outupt.

Next, the program attempts to establish a correspondence between the USNO stars and the stars detected on the image. It does so using a simple (but apparently original) algorithm. The algorithm may not be as elegant as some, but it uses two important pieces of prior information we have here which aren't usually assumed in the general problem, namely (1) we know the image scale in arcsec per pixel quite accurately (from the imagename.hsel file), and (2) we have some idea of the brightnesses of the stars. Thus the algorithm starts with the brightest USNO star, and tentatively assumes that this is the same as the brightest star in the image. It then takes the next brightest USNO star and assumes it's the next brightest star in the image. With two tentative matches, it computes an implied image scale. If this matches the real image scale within half a percent, it starts looking for still more matches. If a third match is found which preserves the image scale with the first two matches, it computes a transformation and starts looking for still more matches. In order to be added to the match list at this point, a star has to fall within 1 arcsec of where it should be. If four or more matches are found, it assumes it's finished and mops up by looking for all further matches. Until it finds four matches, it keeps working its way down the star list, trying different correspondences. This algorithm as coded has something like 6 nested loops and runs formally as N⁶, but in practice it's much faster because it uses the prior information about magnitudes to make good first guesses, and the prior information about the image scale to stay out of the innermost loops until there's a good prospect of success. In most cases a match is found almost immediately. The number of stars in each list is limited to 39, and running down all the possibilities when there are no real matches takes maybe 20 sec on a fast Linux box.

Note that it is quite possible for the algorithm to find four spurious matches and declare victory. This is an insidious error. In single-frame mode, you should examine the output critically to satisfy yourself that this hasn't happened. In batch mode the checkwcs command can be set up to apply several stringent tests to the image set. Read on.

As noted above, once the match is found among the brightest stars, an iterative search is made for correspondences between all the detected stars and all the USNO stars. There's an interesting technical wrinkle here: in the first pass, the parity of the image is unknown, and a 6-constant plate model is used which does not specify the parity of the image (or the perpendicularity of the axes). By comparing the signs of some of the fit parameters, it's possible to deduce the image parity, from this fit. For subsequent fits a four-constant model is used, which requires the image parity to be known a priori. The four-constant model in effect fits a scale factor, a pure rotation, and a zero-point shift - that is, we assume that the CCD has square pixels and perpendicular axes. It's necessary to use this less forgiving fit (once the parity is established) because the situation sometimes arises where all the stars matched are nearly along a diagonal of the picture; because the 6-constant model does not enforce perpendicularity of the axes, the fit then goes bad in the empty corners. It's theoretically possible to have a the initial fit give the wrong parity if the stars are nearly co-linear, but it's more common for the automatch to be completely spurious.

Note that the USNO RA and dec are transformed to tangent-plane coordinates h and x before any fitting is done. Thus there's no worry about distortions near to the pole, etc. [However, the USNO catalog lookup software fails very close to the actual pole, where > 15 min of time is covered by a single picture. Rarely a problem.] You can find a nice discussion of this trick in Taff's Computational Spherical Astronomy, or in any textbook of spherical astronomy. The key is that x and h are local proxies for RA and dec.

Finally, a word about the compression of the USNO used here. At first I'd hoped to use the USNO SA2.0 for this, which is a compressed version of the USNO A2.0. However, the method by which stars were selected for the SA2.0 was evidently to choose stars nearest to specified grid points without regard to magnitude. At low latitude, that means that the USNO SA2.0 looks like `graph paper', and it's nearly impossible to see patterns, either for a human or a machine. Therefore I prepared a compression of the USNO by cutting out everything at d < -45 degrees, and then sifting the remainder by going through 1 square degree at a time and saving the brightest 3000 stars in each square degree. (The r magnitude was used for this; because of the magnitudes are only given to 0.1 mag, I actually cut at a magnitude limit so that 3000 or fewer were kept.) This preserves nearly all the stars at high latitude, but only the brighter objects at low latitude. For the present application, it performs about as well as the full USNO. The compressed version is about 1.15 Gbyte, while the full version is about 6.5 Gbyte, which isn't completely trivial even with today's gigantic disks.

Why matching fails:

I've tested this suite of software extensively, and it works pretty well. But there are lots of ways for it to fail. Here's a checklist:

• Inadequate header info, including no pixel scale, or bad telescope coordinates leading to inadequate overlap with the USNO stars looked up. You'll just have to fix the hsel files somehow. For bad coordinates, you could try the plotusno command to pull up a map of USNO stars in the area, find the pattern, and figure out roughly where you're really pointing. Better to keep your pointing accurate in the first place... If you've got the wrong pixel scale, the CCD software wasn't loaded right. Determine the correct image scale manually (see end of this document), use hedit to fix all the image headers, then try again.

• A bad catalog file (i.e., stars were not properly found in the image). You can see if this is the case by displaying the image from within IRAF, and typing tvmark 1 imagename.cat to mark the stars. If the `stars' are cosmic rays, or don't include the bright stars, or whatever, you may be able to deduce what went wrong. See the cat2seeing discussion for a list of possible pitfalls. It's often fairly straightforward to re-run SExtractor and hand-filter the output to get a better file of stars; see the end of this document for details. Another possibility is to use daofind; this is less convenient for automated operation because it requires an initial guess at the FWHM, but if you're fixing things `by hand' it's reasonably easy.

• Not enough stars in the image. Small-field, short exposure images in insensitive passbands (e.g., U) frequently fail. Also, the USNO has actual holes in it around bright stars, near the Galactic center, and in the middle of bright galaxies and globular clusters - anyplace a Schmidt plate would saturate. There must be four stars matched for the program to think it succeeded. If there are only three good stars, and you're determined, you can use handmatchusno (see end of document) to get a semi-automated match. If there are only two stars, you're out of luck.

• A bad match (Incorrect correspondence between images detected in the picture and the USNO). Although the match algorithm is fairly robust, it can get a spurious match from time to time as noted earlier. The algorithm is not (yet) smart enough to recognize bad matches, `back out' of them, and try again. If you turn on rotation checking in checkwcs later (highly recommended!), nearly all the bad matches will have an implied image rotation out-of-bounds and be caught. But nothing is absolutely foolproof. Again, handmatchusno can be used to patch bad matches.

If at least four matches are found, the matchusno program writes a file called imagename.radec.match.wcs. Each line of this file contains information for a single matched-up star. There's a running number, followed by the x and y coordinates from the picture, followed by the RA and dec from the USNO in `colon-ized' format. This is the main product of the match procedure, and it will be used to feed the IRAF WCS coordinate loader in the next step.

Loading the wcs headers.

In single-picture mode, the matchusno program writes a file wcs_script1 which loads the WCS into the image headers. The safeguard against bad solutions is the user, who must (a) approve the solution before it's written to the script and (b) proceed to run the script. This should be adequately secure.

For batch processing, checkwcs is an intelligent script-writer for this step. It reads the list of filenames; you can also optionally supply a number giving a rotation angle tolerance in degrees, for example

    checkwcs imagelist 1.

Once checkwcs has digested the information for all the images, it makes the following two safety checks: (1) It polls the fits to find the dominant parity. If your data are all taken with the same CCD and telescope, they should always have the same parity! Any image with the opposite parity is rejected. This should kill half the spurious matches. (2) If you've specified a rotation tolerance, it computes the median rotation angle (using an algorithm which is immune to `wrap' problems at zero angle), and rejects any image which is more than the rotation tolerance away from this. This should kill nearly all remaining spurious matches, and occasionally point out cases in which the reference stars are poorly distributed. However, if you adjusted the image rotator partway through the set, this will incorrectly reject fits; in that case you can only apply this test if you parse your data into batches with consistent rotator angle. Images which fail the parity or rotation test are appeneded to the unmatched_pix file.

If very few images are included in the image list, the statistical tests aren't done - in that case a little summary is printed out for a sanity check. If an image comes out at a crazy angle or something, better not use that fit.

checkwcs writes a cl script to the file wcs_script which will do the following for each image with a good solution: (1) Load the wcs information into the image header using the iraf ccmap task. The geometry is specified as rscale, i.e., the same kind of fit used by these programs (parity, rotation, scale, and shift). (2) Use the iraf hedit task to replace the RA and DEC with the computed coordinates of the exact image center. Because (as noted earlier) the WCS procedure commandeers the EQUINOX keyword, this preserves the RA, DEC, and EQUINOX keyword set as a consistent (and now very accurate) description of the pointing.

[Note that checkwcs can also work for a single image, in which case you type checkwsc imagename s. The second argument s tells the program to treat the first argument as the name of an image rather than a list of image names, and suppresses the statistical checks. It looks for the imagename.radec.match.wcs file, fits it, presents the fit to the user, and asks for approval, just as in the last part of matchusno in single-picture mode.]

Checking the results.

checkwcs spews out some hopefully informative output. You can go through the unmatched_pix file to see what may have gone wrong. checkwcs also produces a more detailed file, wcsaudit, which has capsule summaries of each fit (including the number of stars used, the rotation angle, the parity, and the rms). Because both the unmatched_pix file and the wcsaudit file are appended to, rather than created anew each time (as is the wcs_script), the output is time-stamped. That way, if you rerun checkwcs, you can figure out which run of the program started where. Important pictures which didn't match may take some hand work - see the end of this document.

Even Gorier Details

If you want to recompile the code, please make your own copy! If you're contemplating this, I assume you have some idea what you're doing. The codes are in /usr/local/pkg/thor. They link with various things; the makexxx files are not traditional `makefiles' but just executables which compile and link everything which is needed. Utility packages skycalc and numrec are called for, but you shouldn't have to mess with them.

Surgery on the guts of the code is likely to be bloody and painful. There are great swaths adapted from other applications in which the variable names are obscure, etc.; the `double transformation' of going from RA and dec to eta and xi, and then to X and Y, makes for a rather messy business. If you run into some problem I may be able to help (and I'd like to make the code useful to as many as possible); john.thorstensen@dartmouth.edu.

If you want to install the procedures on a machine at home, you'll need the USNO (or better, my reduction of it), PGPLOT for the handmatchusno stuff, and (of course) my various source-code stuff. It shouldn't be hard to do if you're a guru.

On the less esoteric side, you can alter the action of SExtractor by copying the

    /usr/local/pkg/thor/sexdef/starfind.sex

Some uses for the WCS

The WCS stuff is highly integrated into a number of packages. This just scratches the surface.

If you display an image in DS9 (see http://hea-www.harvard.edu) NOT from within IRAF (just display the native FITS file), DS9 knows about the WCS and will show you the RA and Dec of the cursor as you move around. The zoom/align wcs option will put true north at the top. It's also possible to overlay the RA/dec coordinate grid on the image. This is all very cool.

If you display an image from within IRAF on ximtool or DS9, you can run imexamine, and then from within imexamine you type the colon commands

   : set wcs world
   : set xformat %15.3H
   : set yformat %14.2d

Either of these options can be used for a reality check - just graph the USNO stars (perhaps with handmatchusno, display the image with the wcs, and see if the stars' RAs and decs correspond correctly.

The stsdas package for IRAF has some WCS-savvy tools. The imgtools commands under toolbox include rd2xy and xy2rd, wihch do quick conversions (as their names suggest).

The generic IRAF image package, under immatch, has tools wregister and sregister which can be used to align a set of images onto a common grid using the WCS; they can then be mosaicked or averaged. This is a stupendously inefficient way to mosaic images (as far as computer resources go), but it does work. These tasks can be used to register dithered images for stacking; this is a lot easier than laboriously compiling lists of XY shifts for the images you wish to stack.

The USNO A2.0 at MDM

As noted earlier, there's a customized 1.2 Gbyte version of the USNO A2.0 on the Linux boxes at MDM. See the earlier text for details of how this subset was selected. There are a couple of simple utilities to use it:

handmatchusno is a power tool to grab and plot a piece of the USNO and then establish by hand a match between image coordinates and the USNO. It basically does what matchusno does, but with human input to establish the coordinate solution; it's designed to quickly clean up cases where matchusno failed. handmatchusno can also read a pre-existing match file and plot out which stars are used, which is a very useful quick diagnostic. A full description is given later. The underlying graphics package is pgplot

plotusnodisk - reads and plots out a piece of the USNO. The cursor is presented, it lets you find the particulars of the star nearest to the cursor or the coords of the cursor itself. Note that you can direct the graphics to postscript for printing. The underlying graphics package is pgplot. plotusnodisk is operationally a subset of handmatchusno but does not use the hsel information, so may be more general; it also doesn't risk messing with match information.

getusno - Simply reads and dumps out a listing of the USNO.

plotusno - Reads a file from getusno and plots it. Alternatively, you can dump a segment of the full USNO A2.0 from the USNO's website or from Chris Koehn's website at Lowell:

(http://ftp.nofs.navy.mil/data/FchPix/cfra2.html

(http://asteroid.lowell.edu/cgi-bin/koehn/webnet

and plot from them, to fill in all those teeny little missing stars at low latitude.

Note that the USNO has a hole around the Galactic center, and in a few other parts of the sky which was burned out on the original POSS-I plates. Also note that the plates were taken in the 1950s and a (very) few objects have moved significantly since then. Plates north of -15 were POSS-I and quasi-simultaneous; plates south of there are UK Schmidt, are more recent, and were not simultaneous, so the variability can create spurious colors. See the USNO website for greater detail.

There'a hook in handmatchusno which allows you to substitute a faked USNO file for the USNO on the disk. This is helpful in those parts of the sky where the USNO is burned out or otherwise useless. You need to put your star list into the form hh mm ss dd mm ss mag mag. The GSC 1.2 can be useful, but needs to be put into this format. Downloads of the full USNO from the Lowell website are already in this form. See the end of this document for more detail.

Thorstensen's astrometry programs

• handmatchusno is a power tool designed to generate a WCS-friendly input file quickly and easily for the pictures which failed. It's describe in the next section, Fixing Things by Hand.

A couple of simpler programs are as follows:

• astrom - You use this by typing astrom infile, where infile is the name of an input file (the while .radec.match.wcs name must be specified). The program inhales this file, generates a plate solution, and spits out results. The input file format is number X Y RA dec, where number is just an integer label, X and Y are coordinates on the picture, and the RA and dec are the celestial coordinates in sexigesimal format (i.e., hh mm ss dd mm ss). The fields in the celestial coords can be set off either by colons or blanks. If you add objects of known ra and dec to the list with lines of the form number RA dec, it will spit out where they will be in X and Y on the picture. Also, if you put on lines with number X Y, it will spit out the ra and dec of the XY coordinates. So the transformation can be applied either way. The program figures out which is which (i.e., match points for fitting, XY ® ra/dec, or ra/dec ® XY) by noting how many items are on the line, so lines can be in any order. The fit is to the tangent-plane coordinates (good near the poles), and uses a 6-constant plate model for the first fit and a 4-constant (rotation, scale, zero-point) fit for the second pass, after the parity has been determined by the first fit.

• astrompg is the same program as astrom, more or less, but with a graphical display of the fit residuals. I think there's some hook where if you give a negative star number, that residual is plotted but not used in the fit. This is nice when you have a big proper motion.

Fixing Things By Hand...

Sometimes things just don't go right and you have to do something manually. Sometimes it's a matter of fixing up the catalog or hsel file and re-running the automatic matchusno by hand. Other times you have to intervene more closely - for that, I've come up with a power tool for this called handmatchusno. First, a few preliminaries.

I's good to know about the IRAF tvmark command - with this, you can mark on the display the stars which are listed in a *.cat file. To invoke it, pop up an ximtool or ds9 window, and from within IRAF display the image with

 
    display imagename 1

    tvmark 1 imagename.cat

Making your own catalog file. Sometimes if you display and tvmark and image, it'll be obvious that the cat2seeing selection has gone haywire. If you can come up with a decent star catalog, you can try matchusno (the automatic version) again. One easy way to fix this problem might be to adopt an alternate program, like daofind, which requires more input from the user than SExtractor but which works reasonably well, though it will have trouble with saturated stars and requires some parameter tuning for best results. Just rename the output imagename.cat so matchusno will see it.

Another strategy is to filter the output from SExtractor manually. [Note to self: It would be good to write a cursor-driven program to do this interactively.] Start by cutting the line from the matchsript which does the sex command for that image, and re-running it; cut 'n paste will do. That will re-generate the test.cat output file, which you can tvmark again. You can sort the test.cat file using the Linux native sort command. If you wanted (say) to sort in order of the numerical value in the 3rd column, and put the result in a file called temp1, you'd type

    sort -n -k 3 test.cat > temp1

When you're happy, override the bad auto-generated catalog by moving your hand-crafted star catalog file to the be the official *.cat file wich matchusno will look for; for example

    mv temp4 imagename.cat

    !matchusno myimage s

If you're reasonably sure of the fit, apply the wcs to the image by typing

    cl < wcs_script1

Interactive Fitting with handmatchusno

The most stubborn problems occur when the catalog file is good, and the header information is good, but the program somehow gets stuck in an incorrect match. Or perhaps nothing matches because the image scale in the header is wrong, and you don't know what the right value is, so you can't just edit it into the headers. Then you have to figure out the match yourself. The semi-automated program handmatchusno makes this less awkward.

The code will work on its own, but if the .cat and hsel files are available and good, it can make use of them. If the automatic match has failed, these files should exist, so handmatchusno can come along and clean up. But even if they don't exist, it will work with more effort on your part.

Start by displaying your image from IRAF, as described earlier. If a *.cat file exists, tvmark the display for reference, and to verify that it's good. If there is no *.cat file, or if it's bad, you may want to make your own as discussed earlier, but it's not absolutely needed.

Invoke the program from an xterminal outside IRAF by typing

    handmatchusno imagename [seed-imagename] [star-file]

Operation is as follows. The program looks for a *.hsel file, and if it's found it shows you the contents and asks if they're OK. If you say y things proceed; if n (or if it doesn't find the hsel file in the first place), it prompts you for the needed information. The information will be used here only to figure out which piece of the USNO to grab, so you should be reasonably accurate but need not be super-precise. It then asks you if a usable *.cat file exists for this image; answer y or n.

The program looks up the USNO stars in the area (or goes to the star file optionally specified on the command line) and plots them to a graphics window. You're given a cursor, and further commands are cursor-driven:

m

- (The most important command): Match the star nearest the cursor with an XY from the image. After marking the star with m, you're prompted (in the regular xterminal, not graphics window) for the XY coordinates of this USNO star in the image. Get this from the displayed image, possibly with IRAF's imexamine. If a catalog file is in use, the code uses the nearest object in the catalog list, provided it's within 25 pixels of the coordinates you give. In this case you needn't get the XY coordinates exactly right. If you've furthermore specified a seed image, the program shifts the zero-point of that solution using your single match, searches for more matches, and derives a fit if enough matches are found. If you're not using a catalog file of detected images, the match is simply recorded, so you need to provide a fairly precise XY, perhaps from imexam.

o

- Overplot little boxes around all the matched stars, with XY coordinates labeling them. Useful for checking solutions.

b

- Plots a box showing the boundaries of the picture; only works if there are enough stars for a valid solution.

f

- Fit the matches we have so far. At least three matches are needed to make a fit. If your matches are valid, the RMS error of the fit will generally be less than 1 arcsec. If you don't have a good catalog you'll want to enter at least four matches by hand with m before trying a fit, in order to get some redundancy.

r

- `refine' - Using the fit we have, go looking for more matches. There must be a catalog file for this to work.

s

- `sigma-clip' - delete matches whose residuals are more than a specified number of times larger than the rms residual, and re-fit.

e

- `edit' - list the matches and the fit errors in RA and dec (arcsec), and optionally delete one of the matches (sorry, this is primitive and only lets you kill one match at a time. See the s option above.)

g

- Gets an existing match file for this field (imagename.radec.match.wcs suffix), and computes the solution for it; you can then use o to overplot the matches. This is useful for checking the reliability of a previously-existing solution. If you find that the match file has some flaw (e.g., a badly-centered saturated star has pulled the fit), you may be able to fix it using e, m, and r, which manipulate the contents of the match file.

p

- (re)Plot the stars but not the little boxes (fresh plot).

q

- Quit, writing results to a radec.match.wcs file and generating a wcs_script1 script to load the coordinates.

x

- Exit, without writing out results.

    cl < wcs_script1

Using handmatchusno with a star file.

In normal operation handmatchusno looks up the catalog stars from the USNO on disk. As explained earlier, you can override this and use your own file of catalog stars by giving a third command line argument (the first and second being the filename roots of the program frame and a seed calib frame), which is the name of your file of catalog stars. The format expected is hh mm ss dd mm ss mag1 mag2; anything can be floating point, and the format is otherwise flexible. mag1 is normally red and mag2 is blue, but this isn't critical. Using an alternate star list works even if you don't have a seed calibration; just put in xxx or something as a placeholder for the second argument. The program will try and fail to look up a seed solution, but will then proceed normally.

There are several sources for star files. If the truncated USNO is too sparse (which happens rarely) you can download a properly-formatted excerpt of the full USNO from

http://asteroid.lowell.edu/cgi-bin/koehn/webnet

http://www-gsss.stsci.edu/gsc/gsc2/GSC2home.htm

   reformgsc2 filename

The USNO A2.0 and GSC catalogs have somewhat limited accuracy. As an experiment I obtained the UCAC-1 catalog, which is the USNO's very precise CCD astrometric catalog covering most of the southern hemisphere. I haven't installed it at MDM, but am describing my experience here to make things smoother for others who may wish to use it. The UCAC-1 positions should be at least three times better than the USNO's, and it does better in the Galactic center, but not as well as GSC2.2. In one low-latitude test I found the UCAC-1 to be missing so many stars (rejected because of crowding) that I couldn't easily match the patterns, but since I had a solution based on the GSC it wasn't hard to find establish the single matching star needed to jump to the solution. As expected, the UCAC-1 solution had very low residuals (like 0.1 arcsec) and matched lots of stars in the image (typically 100 in an 8x8 arcmin field near the Galactic center). I haven't set up software for automatically reading the UCAC-1 in handmatchusno - instead I used the software supplied with the catalog to manually read in the stars in the region, and another little piece of code to transform the formats to my expected USNO list format. This is too much work for routine use, but is fine for the occasional patch. I'm happy to supply little bits of the catalog by request; if you need it on a routine basis you'll need to get the CD-ROM from USNO and set it up yourself (the binary file needs 600 Mbyte and doesn't truncate gracefully without some work on your part).