Synopsis
Create a frame transformation to minimize the aspect difference between data from the same sky region.
Syntax
wcs_match infile refsrcfile outfile [wcsfile] [logfile] [method] [radius] [select] [residlim] [residtype] [residfac] [multimatch] [clobber] [verbose]
Description
wcs_match compares two sets of source lists from the same sky region. If three or more sources are found to be a close match in position, wcs_match will calculate a transformation matrix which, when used with wcs_update, can be used to align the input (infile) set of source positions with the reference (refsrcfile) set of source positions. If only one source is found to be a match, a simplified, translationonly transformation matrix can be computed.
The position differences are minimized in a leastsquares sense, minimizing the sum of the squared residuals. More details of the method used can be found in Appendix B.2 of Practical Handbook on Image Processing for Scientific Applications by Bernd Jahne.
Options allow the user to specify how close sources from each list must be to each other to be considered as possible matches (radius) and to limit the acceptable error between reference source positions and updated input source positions (residlim, residtype and residfac). The calculated six element transformation matrix is output to the transform file (outfile).
WCS parameters specifying a tangent point must be provided for the transform calculations. This can be done by specifying either a FITS image or table file containing a WCS in the wcsfile input parameter. Alternately, if a WCS is already present in either the input or reference source lists, those values will be used and wcsfile can be left blank. When a wcsfile is specified, it should contain a WCS that identifies a tangent point that is relatively close to the sources to be matched.
Transform Calculation Controls
wcs_match calculates a transformation matrix which minimizes (in a least squares sense) the error between tangent plane projections of the reference sources and the transformed input sources. Using default method=rst (rotate, scale translate), the transformation matrix specifies a 2D translation, a rotation about the tangent point, and a scale factor. Applying this matrix to the input sources results in the best attempt at aligning the reference and transformed input source positions. A simpler, translationonly transformation matrix is available by changing the method parameter.
A transformation matrix is initially calculated using data from all matched input and reference source positions (source pairs). The transformation matrix is then applied to the input source positions and the updated positions are compared to the positions of the reference sources. The differences between these source positions (the residuals) are used to determine whether the transformation matrix is sufficient, or whether certain source pairs should be excluded (e.g., mismatched sources) and the transformation calculation repeated. Running wcs_match with verbose=1 shows details of the residuals during successive calculations of the transformation matrix.
If residlim is greater than 0, the transform calculation stops when the largest residual is smaller than residlim. If the largest residual is greater than residlim, the source pair with this residual is excluded from the matched set of sources and the transformation matrix is recalculated. This continues iteratively until all residuals are less than the specified value of residlim.
If residfac is greater than 0, the residualtosource pair position error ratio is examined. This ratio is determined by dividing a residual by the sum in quadrature of its reference and input source position errors. A comparison is made to residfac and if residfac is exceeded, the source pair with the largest ratio is excluded from the matched set of sources and the transformation matrix is recalculated. This continues iteratively until residfac is no longer exceeded. Two different comparisons can be made with residfac. If residtype = 0, the transform calculation stops when all ratios are less than residfac. If residtype = 1, the calculation stops when the average of all ratios is less than residfac.
Sources are matched if they are within the specified radius. By default if there are 2 or more sources matching the same input source or reference source, then that source is excluded from the calculation. That is only unique, unambigious matches are considered. New in CIAO 4.14, users can set the multimatch parameter to yes, which will allow sources matching multiple input or reference sources to be included in the calculations.
Examples
Example 1
wcs_match src1.fits src2.fits xfm.fits wcsfile=image.fits radius=2 residlim=1 verbose=1
Match sources from src1.fits and src2.fits that are within 2 arcseconds. Eliminate any source pairs with a position difference greater than 1 arcsecond from the transform calculation. Output the elements of the inputtoreference source position transformation matrix in xfm.fits. Display the intermediate and final position errors for all matched sources (verbose=1).
Example 2
wcs_match src1.fits src2.fits xfm.fits method=trans verbose=0
Use default values for matching sources and limiting source pairs used in transform calculation. Use a WCS for the transform calculation from either the input source file (infile) or the reference source file (refsrcfile). Create a transformation matrix which uses only translation, not scale or rotation, to match source positions (method=trans). Do not display any source pair position error information (verbose=0).
Example 3
wcs_match src1.fits src2.fits xfm.fits wcsfile=image.fits residlim=0 residtype=1 residfac=1.5 verbose=1
Do not use the absolute value of the residual to limit source pairs to include in the transform calculation (residlim=0), but do use the residualtoposition pair position error ratio to limit source pairs in the transform. Whenever the average of all ratios (residtype=1) is greater than 1.5, remove the source pair with the largest ratio from the calculation and repeat. Display the intermediate and final position errors for all matched sources (verbose=1).
Parameters
name  type  ftype  def  min  max  units  reqd 

infile  file  input  yes  
refsrcfile  file  input  yes  
outfile  file  output  yes  
wcsfile  file  input  no  
logfile  file  output  STDOUT  no  
method  string  rst  no  
radius  real  12  arcseconds  no  
select  string  auto  no  
residlim  real  5  arcseconds  no  
residtype  integer  0  no  
residfac  real  2.0  no  
multimatch  boolean  no  no  
clobber  boolean  no  no  
verbose  integer  0  0  5  no 
Detailed Parameter Descriptions
Parameter=infile (file required filetype=input)
File with input sources.
File with input source positions. Must contain columns RA and Dec., either internally or through use of datamodel column filter. RA and Dec. columns must be in decimaldegree (i.e., not sexagesimal) format. At least 3 input sources must be found to be a match for 3 reference sources. residfac parameter will only be useful if columns RA_err and Dec_err are present in infile and refsrcfile.
Parameter=refsrcfile (file required filetype=input)
File of reference sources
File of reference source positions. Must contain columns RA and Dec., either internally or through use of datamodel column filter. RA and Dec. columns must be in decimaldegree (i.e., not sexagesimal) format. At least 3 input sources must be found to be a match for 3 reference sources. residfac parameter will only be useful if columns RA_err and Dec_err are present in infile and refsrcfile.
This file can be in any format supported by the CXC Datamodel, including simple ASCII tables, such as
#RA RA_ERR DEC DEC_ERR 246.8147625543 1.6359623089102e05 24.55049151732 1.7897900878694e05 246.8861311412 1.3138809151769e05 24.55677209605 9.1986712043024e06 246.8794236861 1.4607012701617e05 24.56762593616 1.0942152091076e05 246.8747572911 2.1446620706911e05 24.56022597575 1.4661890684664e05 ...
This may be especially convenient when matching a Chandra dataset to data from another waveband (optical, IR) or external catalog.
Parameter=outfile (file required filetype=output)
File containing transform matrix data.
This file contains the six elements of the transform matrix and WCS coordinate information used in calculating the transform. It can be used by wcs_update to update FITS images, tables, or asol files.
Parameter=wcsfile (file not required filetype=input)
World Coordinate System (WCS) file.
The WCS file must contain WCS coordinate transform data with a tangent point near the sources of the infile or refsrcfile. It can be either a FITS image file with a "EQPOS" WCS, or a FITS table file with "EQPOS" or "EQSRC" WCS. A WCS file must be specified unless a WCS is present in either the input or reference source file. If a wcsfile is provided and a WCS is present in either source file, the WCS from the wcsfile will be used.
Parameter=logfile (file not required filetype=output default=STDOUT)
Debug log file.
Allowable values are either a filename  to send the output to a given file  or one of stdout or STDOUT, which sends the information to the standard output (normally the screen).
Parameter=method (string not required default=rst units=)
Method used to generate transformation matrix
The transformation matrix used to align input source positions with reference source positions can be caculated in two ways. If method='rst', the transformation matrix includes a rotation, scale factor, and translation. To limit the transformation matrix to provide only a translation, method can be set to 'trans'.
Parameter=radius (real not required default=12 units=arcseconds)
Source match radius
Source positions from the two source files can only match each other if they are no farther apart than the value specified as the source match radius. Only a single reference source can be within the source match radius of a input source, and only a single input source can be within the source match radius of a reference source. If multiple sources are found in either case, the source pair is not used to calculate the transform.
Parameter=select (string not required default=auto)
Automatically or manually select which source pairs to use to compute the transform. Valid options are "auto" and "manual".
With select="auto", the tool will automatically determine the best set of matching sources based on the earlier descriptions of the algorithm. With select="manual", the tool will enter an interactive mode where the user can iteratively select which source pairs to include or exclude when computing the transformation.
When method="rst", the tool still requires at least 3 sources; even after manual source selection. However, with method="trans" if all the source pairs are excluded, the output will have a "unit" transform (0 rotation, unit scale, and 0 translation).
Parameter=residlim (real not required default=5 units=arcseconds)
Residual Limit
After the transform is calculated and initially applied to the input source positions, if the largest source pair position error exceeds the value of residlim, then that pair is eliminated from the transform calculation and the calculation is repeated. This process of eliminating the source pair with the largest residual and repeating the calculation continues until all remaining source pair position errors are less than the residlim value.
Set residlim to 0 to disable this feature and prevent source pairs from being removed based on absolute residual value. Residlim takes precedence over residfac if both are nonzero.
Parameter=residtype (integer not required default=0 units=)
Residual ratio limit type
Two different options can be used to eliminate source pairs from the transform calculation when the residual ratio limit factor (residfac) is specified. If residtype is set to 0, residfac is interpreted as an upper limit on the residualtosource pair position error ratio. Any sources with ratios above residfac will be eliminated from the transform calculation. If residtype is set to 1, residfac is interpreted as an upper limit on the average of all ratios. Source pairs with the highest ratios will be incrementally removed from the transform calculation until the average of all ratios drops below residfac.
Parameter=residfac (real not required default=2.0 units=)
Residual ratio limit factor
After the transform is calculated and initially applied to the input source positions, if either the largest residualtosource pair position error ratio exceeds the value of residfac (when residtype=0), or the average of all ratios exceeds the value of residfac (when residtype = 1) then the source pair with the largest ratio is eliminated from the transform calculation and the calculation is repeated. This process of eliminating a source pair and repeating the calculation continues until the appropriate ratio test is met.
Set residfac to 0 to disable this feature and prevent source pairs from being removed based on the residualtosource pair position error ratio. Residlim takes precedence over residfac if both are nonzero.
Parameter=multimatch (boolean not required default=no)
Allow multiple input or reference sources to match?
By default only unique source matches are considered when computing the transform matrix. By setting multimatch=yes, users can allow sources that match multiple input or reference sources to be included in the transformation computation.
Parameter=clobber (boolean not required default=no)
Overwrite existing output dataset with same name?
Parameter=verbose (integer not required default=0 min=0 max=5)
Level of debug detail.
Increasing amounts of debug information is printed to "logfile" as the value of verbose is increased from 0 to 5. Setting verbose=1 shows details of the source pair errors after applying the transform to the input source positions.
Interactive mode with select='manual'
In CIAO 4.13, users can now manually, interactively select which source pairs are to be used when computing the transform matrix. This is done by setting the parameter select='manual'.
In this mode, the tool will compute the offsets of the matching sources and will then prompt the user the supply a list of source pairs to include or exclude. The list is used to invert the current state of a source pair: so if the source pair is included it becomes excluded; and if it is excluded then it becomes included. The residuals and matrix are updated after each update and the user can repeat the process. The iterations end when the user inputs "1" (no quotes).
For example consider this example
% wcs_match 18898.dat 20056.dat out=mytransform.fits wcs=acisf20056_repro_evt2.fits \ select=manual clob+ method=trans Source Residuals  Src Ref# Dup# Ref RA Ref Dec. Prior Resid Transfm Resid Resid Incl Indx (deg.) (deg.) RSS (x,y) RSS (x,y) Ratio (arcsec) (arcsec) 0 0 0 186.71376 33.27030 0.57 ( 0.27, 0.51) 0.00 (0.00,0.00) 0.00 Y Source Residuals, before/after transform (arcsec), and percentage improvement: Average Residuals: 0.574413 0.000000 100.00% Maximum Residuals: 0.574413 0.000000 100.00% RMS Residuals: 0.406171 0.000000 100.00% Source Residual Ratios, before/after transform, and percentage improvement: Average Residual Ratios: 0.583753 0.000000 100.00% Maximum Residual Ratios: 0.583753 0.000000 100.00% RMS Ratios: 0.412775 0.000000 100.00% Please enter the Src Indx number of a source to delete (or add back), or a comma separated list of Src Indx numbers to delete (or add back). Enter a 1 to finish: 0
We see in this example that there is 1 matching source pair. Taking a look at the "Source Residuals" table, we see that the first column is the source index number, Indx, is 0, and that the last column, Incl, indicates that yes, "Y", the source pair is being included.
By entering '0' when prompted (no quotes), we are asking to now exclude source pair 0. The updated Source Residuals table is then shown along with the updated residuals. Since we are using method=trans, and 0 sources pairs are left, the tool is using a "unit" transform (no rotation, unit scale, no translation).
Source Residuals  Src Ref# Dup# Ref RA Ref Dec. Prior Resid Transfm Resid Resid Incl Indx (deg.) (deg.) RSS (x,y) RSS (x,y) Ratio (arcsec) (arcsec) 0 0 0 186.71376 33.27030 0.57 ( 0.27, 0.51) 0.57 ( 0.27, 0.51) 0.58 N Source Residuals, before/after transform (arcsec), and percentage improvement: Average Residuals: 0.574413 0.574413 0.00% Maximum Residuals: 0.574413 0.574413 0.00% RMS Residuals: 0.406171 0.406171 0.00% Source Residual Ratios, before/after transform, and percentage improvement: Average Residual Ratios: 0.583753 0.583753 0.00% Maximum Residual Ratios: 0.583753 0.583753 0.00% RMS Ratios: 0.412775 0.412775 0.00% Please enter the Src Indx number of a source to delete (or add back), or a comma separated list of Src Indx numbers to delete (or add back). Enter a 1 to finish: 1
To end the interation users type "1" (no quotes).
As of CIAO 4.16, users can now specify a range of values or input a stack (external file) . A range example:
0:3
will flip the status of Indx values 0, 1, 2, and 3. An example using stacks:
@my_file.lis
where the file named my_file.lis contains the values 1 and 2, would flip the status of those two Indx values.
Changes in the CIAO 4.16 release
Improved logic for deleting poor matches
The tool will remove sources from one list which match multiple sources in the other list. If the last matching pair hit this condition then in CIAO 4.15 and earlier the tool would have deleted both this match and the preceeding (valid) match. The valid match is nolonger deleted.
Manual source selection
Sources can now be specified as a range  with the format "low:high"  or a stack file. If a list contains repeated sources then the query will be repeated.
Changes in CIAO 4.16

If the last matching pair was a poor match, the code mistakenly deleted that pair and the prior, valid matched pair. The logic has been updated to prevent the valid pair from being deleted.

Fixed a bug when using manual matching and duplicate values were input.

In manual matching mode, users can now specify values in an external file using the @file_name syntax as well as supporting ranges in the form a:b.
Bugs
Caveats
 Using method=rst (default) may change the pixel size/scale which introduces problems for some downstream analysis.

When running wcs_match with the default parameter setting: method=rst (rotate, scale, translate), it may compute an astorometric solution that includes a change to the pixel size (scale). Typically this is then applied to the Level2 event file. The updated event list may then have a nonstandard pixel size which is propagated into images created by binning the event file.
However, some data products such as exposure maps, are always created using the standard pixel sizes. This mismatch in pixel size can then trigger errors or unexpected results when trying to use the two data products together.
There are two types of errors. (1) Using the same grid to bin the event file and create the exposure map may result in different image axis lengths; typically they will be off by 1 pixel. (2) Attempting to use the exposure map as a Datamodel mask() filter, will result in an error since it requires the mask and the image being filtered to have the same pixel size; this type of filtering is done with merge_obs when creating a combined PSF map.

Workaround:
Users should only run wcs_match using method=trans which forces the astrometric solution to only include linear translations.
 Using method=rst (default) may result in a nonphysical astrometric solution.

When running wcs_match with the default parameter setting: method=rst (rotate, scale, translate), it may compute an astorometric solution that is nonphysical. Examples of nonphysical solutions include large translation offsets which are compensated for with large rotation offsets in the opposite direction. The Correcting Absolute Astrometry thread has example of how to interpret wcs_match's output transform file. Given the accuracy of the Chandra pointing, offsets larger than a few pixels are highly suspicious. In addition, since the aspect camera field of view is much larger than the instruments' field of view, rotation adjustments larger than a fraction of a degree are unlikely physical.
This is caused by having relatively few matching sources and/or all the matching sources are clustered together on one side of the tangent point. Mathematically, to use method=rst requires at least 3 source pairs, however, in practice users should have dozens of matching sources dispersed throughout the entire field of view. With few sources, the underconstrained "rst" solution may be mathematically optimal, but physically impossible.

Workaround:
Users should only run wcs_match using method=trans which forces the astrometric solution to only include linear translations.
See Also
 tools::coordinates
 reproject_aspect, reproject_events, sso_freeze, wcs_update