wcs_match

Synopsis

Create a frame transformation to minimize the aspect difference between data from the same sky region.

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, translation-only transformation matrix can be computed.

The position differences are minimized in a least-squares 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 2-D 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, translation-only 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 residual-to-source 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 input-to-reference 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 residual-to-position 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

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 decimal-degree (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 decimal-degree (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.6359623089102e-05 -24.55049151732 1.7897900878694e-05
246.8861311412 1.3138809151769e-05 -24.55677209605 9.1986712043024e-06
246.8794236861 1.4607012701617e-05 -24.56762593616 1.0942152091076e-05
246.8747572911 2.1446620706911e-05 -24.56022597575 1.4661890684664e-05
...

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 residual-to-source 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 residual-to-source 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 residual-to-source 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 no-longer 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 down-stream 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 non-standard 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 non-physical astrometric solution.

When running wcs_match with the default parameter setting: method=rst (rotate, scale, translate), it may compute an astorometric solution that is non-physical. Examples of non-physical 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 under-constrained "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