Last modified: 22 Feb 2022

URL: https://cxc.cfa.harvard.edu/ciao/threads/reproject_aspect/

Correcting Absolute Astrometry

CIAO 4.16 Science Threads


Overview

Synopsis:

Chandra's absolute pointing accuracy is generally better than 0.4 arcsec. This can be improved in an absolute sense by cross-matching Chandra sources with other higher precision catalogs or in a relative sense by cross matching one Chandra observation with another observation. The CIAO tools reproject_aspect, wcs_match, and wcs_update can be used to compute the fine astrometic shifts between two source lists and apply the offsets to various Chandra files.

reproject_aspect is actually a script which runs the two tools: wcs_match and wcs_update. For more flexibility these tools can be run individually which is done in this thread.

Purpose:

  • To reduce the errors in the absolute or relative astrometry in a Chandra dataset.

Related Links:

Last Update: 22 Feb 2022 - Review for CIAO 4.14. Updated for Repro5/CALDB 4.9.6


Contents


Get Started

Download the sample data: 2313 (Chandra Deep Field South); 12220 (Chandra Deep Field South)

unix% download_chandra_obsid 2313,12220 

In this thread, we assume that all relevant files are in the same working directory.

The two Chandra observations have very different observation times

unix% dmkeypar acisf02313_repro_evt2.fits ONTIME echo+
132124.93587506
unix% dmkeypar acisf12220_repro_evt2.fits ONTIME echo+
48768.153534293

In this example the longer observations, OBS_ID = 2313, will be used as the reference dataset. The assumption will be that the longer observation will have better statistics and produce better source centroids. The fine astrometric shifts will be applied to the shorter observation, OBS_ID = 12220 to register it to match the longer observation.


Matching Sources

There are two common starting places to perform fine astrometric corrections. Either by cross matching Chandra sources with an external catalog whose positions are expected to be better than the Chandra positions, or by cross matching one Chandra observation with another Chandra observation which provides a relative correction.

This thread will show an example of both.

This thread will use the individual tools wcs_match and wcs_update which make up the reproject_aspect script.

Running wavdetect

To begin with, wavdetect will be used to detect sources in Chandra observation ObsID 12220. Users can review the Running wavdetect thread for a full discussion of the wavdetect parameters and how to create the various input files. A summary of the commands used is shown here:

Create counts image and exposure map
unix% fluximage acisf12220_repro_evt2.fits cdfs12220 bin=1 band=broad

Create PSF map
unix% mkpsfmap cdfs12220_broad_thresh.img cdfs12220_broad_thresh.psfmap energy=2.3 ecf=0.9 mode=h clob+

Run wavdetect
unix% punlearn wavdetect
unix% wavdetect \
  infile=cdfs12220_broad_thresh.img \
  psffile=cdfs12220_broad_thresh.psfmap \
  expfile=cdfs12220_broad_thresh.expmap \
  scales="1 2 4 6 8 12 16 24 32" \
  outfile=cdfs12220_broad.src \
  scell=cdfs12220_broad.cell \
  imagefile=cdfs12220_broad.recon \
  defnbkg=cdfs12220_broad.nbkg \
  interdir=./ mode=h clob+

Briefly: fluximage is used to create the Chandra Source Catalog (CSC) broad band, 0.5 - 7.0keV, image and 2.3 keV monochromatic exposure map. mkpsfmap is used to create the corresponding PSF map; a large ecf=0.9 value is used to help suppress small faint sources. For the purposes of registering sources, only good quality detections with robust centroids should be used. Finally, the wavdetect tool is used to locate sources. Figure 1 is the broad band image showing the wavdetect regions.

[NOTE]
Changes in CIAO 4.11: PSF map file

The fluximage script now creates the PSF map file when the user supplies the psfecf value at the same energy as used to create the exposure map.

unix% pset fluximage psfecf=0.9

Users can also continue to independently run mkpsfmap.

Figure 1: Wavdetect source regions for ObsID 12220

[Thumbnail image: Wavdetect source regions for ObsID 12220]

[Version: full-size]

[Print media version: Wavdetect source regions for ObsID 12220]

Figure 1: Wavdetect source regions for ObsID 12220

Broad band, 0.5 to 7.0 keV, counts image for ObsID 12220. The green ellipses are the wavdetect source detections.

In general users should further scrutinize and filter the wavdetct output source lists (filter on source significance, or net counts); however for brevity this step is omitted from this thread.

Users do not need to specifically run wavdetect. Any source detection tool including celldetect and vtpdetect or any number of other non-CIAO tools may be used.


External Catalog

The first example will show how to match the Chandra observation with an external catalog.

The choice of which catalog is strongly driven by the science being done. Different catalogs in different wavelengths cover different parts of the sky with different sensitivity and different accuracy.

Since ds9 already has an interface to various standard catalogs this thread will use it to retrieve data from the United States Naval Observatory, USNO-A2.0, catalog . The USNO-A2.0 catalog can be queried and displayed from the command line using the -catalog switch

unix% ds9 cdfs12220_broad_thresh.img -catalog ua2

or it can be accessed by going to Analysis → Catalogs → Optical → USNO-A2.0. The catalog is displayed in Figure 2.

Figure 2: ObsID 12220, overlaid with USNO-A2.0 optical catalog.

[Thumbnail image: ObsID 12220, overlaid with USNO-A2.0 optical catalog.]

[Version: full-size]

[Print media version: ObsID 12220, overlaid with USNO-A2.0 optical catalog.]

Figure 2: ObsID 12220, overlaid with USNO-A2.0 optical catalog.

The broad band, 0.5 to 7.0 keV, counts image for ObsID 12220 overlaid with the USNO-A2.0 optical catalog (light blue circles). The position of the wavdetect source locations are shown as red, plus ("+") symbols.

Note: The size of the symbol is arbitrary and not to any scale; it does not represent the uncertainty or error in the catalog positions.

The USNO-A2.0 catalog can now be saved to an ASCII file. In the Catalog window, goto File → Export → Tab-Separated-Values ... and save the file as usno_a2.tsv in the same directory as the other files used in this thread. Inspecting the ASCII catalog file:

unix% head -5 usno_a2.tsv
_RAJ2000	_DEJ2000	USNO-A2.0	RAJ2000	DEJ2000	ACTflag	Mflag	Bmag	Rmag	Epoch
052.77165300	-27.84893100	0600-01389895	052.771653	-27.848931			19.6	17.8	1983.246
052.77216200	-27.81999500	0600-01389915	052.772162	-27.819995			18.9	17.1	1983.246
052.77362000	-27.81208100	0600-01389964	052.773620	-27.812081			20.3	18.1	1983.246
052.77406400	-27.86205300	0600-01389981	052.774064	-27.862053			20.4	18.1	1983.246

shows that it is not one of the standard formats that the DM recognizes so it requires the use of the [opt ] qualifier to skip the header line and jump directly to the columns of data. This can be verified using dmlist

unix% dmlist "usno_a2.tsv[opt skip=1]" cols

--------------------------------------------------------------------------------
Columns for Table Block usno_a2.tsv
--------------------------------------------------------------------------------
 
ColNo  Name                 Unit        Type             Range
   1   col1                              Real8          -Inf:+Inf            
   2   col2                              Real8          -Inf:+Inf            
   3   col3                              String[13]                          
   4   col4                              Real8          -Inf:+Inf            
   5   col5                              Real8          -Inf:+Inf            
   6   col6                              String[4]                           
   7   col7                              Real8          -Inf:+Inf            
   8   col8                              Real8          -Inf:+Inf            

wcs_match takes in these two catalogs (source lists), performs a simple cross match, and determines the transformation parameters to shift input sources to the reference source locations in a way that minimizes the difference. This process is iterative and may exclude some cross matches if by removing them the solution is improved.

wcs_match can either determine a solution that is translation only (shifts in RA and Dec.), method="trans" or it can compute a solution that involves rotation, plate-scale, and translation corrections, method="rst" (default). While more general, method="rst" should be used carefully and only when there are many pairs of matching sources evenly distributed throughout the field. Since this example has few nearby, matching sources it will use method="trans" to determine a simple translation only solution.

The wcs_match parameters are now set as follows

unix% pset wcs_match infile=cdfs12220_broad.src
unix% pset wcs_match refsrcfile="usno_a2.tsv[opt skip=1][cols ra=col1,dec=col2]"
unix% pset wcs_match outfile=out.xform
unix% pset wcs_match wcsfile=cdfs12220_broad_thresh.img 
unix% pset wcs_match method=trans

Here the DM columning renaming syntax, [cols ra=col1,dec=col2] is used to identify the RA and Dec columns in the ASCII file. The wcsfile is set to the image file for ObsID 12220. Users may want to set some of the other hidden parameters that control the cross match algorithm that include the search radius and convergence criteria. After the parameters are set, the tool is executed

unix% wcs_match 
Input file with duplicate srcs (cdfs12220_broad.src): 
Input file with reference srcs (usno_a2.tsv[opt skip=1][cols ra=col1,dec=col2]): 
Transform file (out.xform): 
# wcs_match (CIAO): WARNING: Dup src x_err or ra_err cols not found.  Assuming x_err = 1.
# wcs_match (CIAO): WARNING: Dup src y_err or dec_err cols not found.  Assuming y_err = 1.

The warnings are expected since the reference source list, usno_a2.tsv, does not contain estimates of the position errors. Running the tool with verbose > 0 will provide information about which sources were used in determining the final solution.

The output file, out.xform, looks like

unix% dmlist out.xform cols
 
--------------------------------------------------------------------------------
Columns for Table Block out.xform
--------------------------------------------------------------------------------
 
ColNo  Name                 Unit        Type             Range
   1   a11                  n/a          Real8          -Inf:+Inf            Transform Matrix element A11
   2   a12                  n/a          Real8          -Inf:+Inf            Transform Matrix element A12
   3   a21                  n/a          Real8          -Inf:+Inf            Transform Matrix element A21
   4   a22                  n/a          Real8          -Inf:+Inf            Transform Matrix element A22
   5   t1                   n/a          Real8          -Inf:+Inf            Transform Matrix element T1
   6   t2                   n/a          Real8          -Inf:+Inf            Transform Matrix element T2
   7   ra_ref               deg          Real8          -Inf:+Inf            WCS tangent point RA
   8   dec_ref              deg          Real8          -Inf:+Inf            WCS tangent point Dec.
   9   roll_ref             deg          Real8          -Inf:+Inf            WCS tangent point roll
  10   xpix_ref             pixel        Real8          -Inf:+Inf            WCS tangent point x (pix)
  11   ypix_ref             pixel        Real8          -Inf:+Inf            WCS tangent point y (pix)
  12   x_scale              deg-per-pixel Real8          -Inf:+Inf            WCS scale in x direction
  13   y_scale              deg-per-pixel Real8          -Inf:+Inf            WCS scale in y direction

The first 6 columns are the most interesting. They provide the matrix of coefficients needed to shift the input source positions to match the reference source locations.

\[ \left( \begin{array}{c} X'\\ Y' \end{array} \right) = \left( \begin{array}{ll} a_{11} && a_{12} \\ a_{21} && a_{22} \end{array} \right) \left( \begin{array}{c} X \\ Y \end{array} \right) + \left( \begin{array}{c} t_1 \\ t_2 \end{array} \right) \]

where in terms of a scale factor, S and a rotation angle, θ

\[ \begin{array}{l} a_{11} \approx S*\cos\theta \\ a_{12} \approx S*\sin\theta \\ a_{21} \approx -S*\sin\theta \\ a_{22} \approx S*\cos\theta \end{array} \]

The values are approximately equal since the general solution does allow for skew; however, when applied to the data, just a rotation is applied.

Examining the output values

unix% dmlist out.xform"[cols a11,a12,a21,a22,t1,t2]" data,clean
#  a11                  a12                  a21                  a22                  t1                   t2
                  1.0                    0                    0                  1.0    -0.19426639288812     0.49610366114037

Since method=trans was use, the a coefficients are either 0 or 1.

[IMPORTANT]
Units

It is important to point out that the coefficients are in units of the physical pixels of the wcsfile. This is why it is critical to use the same wcsfile in both wcs_update and wcs_match

This part of the thread is done. The output file out.xform can now be used in the Apply transformation to Chandra Dataset section of this thread.


A Second Chandra Observation

This section can be skipped if matching with an external catalog (above).

The steps for using a second Chandra observation as the reference source list is basically the same as the External Catalog, except the reference source list will just happen to be a second Chandra datasets. This computes the offsets needed to register the two Chandra datasets, which is often the first step prior to merging.

As before, for a full discussion on running wavdetect users should review the Running wavdetect thread. The steps used are shown below to generate a source list for the second Chandra observation, ObsID 2313.

Create counts image and exposure map
unix% fluximage acisf02313_repro_evt2.fits cdfs2313 bin=1 band=broad

Create psfmap
unix% mkpsfmap cdfs2313_broad_thresh.img cdfs2313_broad_thresh.psfmap energy=2.3 ecf=0.9 mode=h clob+

Run wavdetect
unix% punlearn wavdetect
unix% wavdetect \
  infile=cdfs2313_broad_thresh.img \
  psffile=cdfs2313_broad_thresh.psfmap \
  expfile=cdfs2313_broad_thresh.expmap \
  scales="1 2 4 6 8 12 16 24 32" \
  outfile=cdfs2313_broad.src \
  scell=cdfs2313_broad.cell \
  imagefile=cdfs2313_broad.recon \
  defnbkg=cdfs2313_broad.nbkg \
  interdir=./ mode=h clob+

The wavdetect source locations are shown in Figure 3.

Figure 3: CDFS Source Wavdetect Source Positions

[Thumbnail image: CDFS Source Wavdetect Source Positions]

[Version: full-size]

[Print media version: CDFS Source Wavdetect Source Positions]

Figure 3: CDFS Source Wavdetect Source Positions

The broad band counts image for ObsID 12220 with the wavdetect source locations plotted. The red plus, "+", crosshairs are the source detections from this ObsID. The blue "x" are from the longer dataset, ObsID 2313. As expected, there are many more sources detected in the longer dataset.

Note: The symbol size is arbitrary. It does not represent the error or uncertainty in the source positions.

These two source list can now be input to the wcs_match tool to compute the astrometric correction needed to register the shorter observation, ObsID 12220 to the longer observation ObsId 2313.

unix% pset wcs_match infile=cdfs12220_broad.src
unix% pset wcs_match refsrcfile=cdfs2313_broad.src
unix% pset wcs_match outfile=out.xform
unix% pset wcs_match wcsfile=cdfs12220_broad_thresh.img
unix% pset wcs_match method=trans

As discussed above, method=trans is used to restrict the solution to only include a translation. For this specific example there are a large number of sources which could yield reliable rotation and scale components; however, for consistency with the earlier only the translation is used here.

[NOTE]
wcsfile

One common question is what to use for the wcsfile parameter. The answer is it does not matter as long as you use the same wcsfile for both wcs_match and wcs_update. The tools work by making the corrections in a tangent plane. The choice of that tangent plane is arbitrary; it just needs to be consistent between the two tools. It should be close to the location of the being examined but does not specifically need come from either dataset.

In this example, the shorter observation is used for the wcsfile even though it is the one being corrected. The final results are not affected.

Now run the tool

unix% wcs_match
Input file with duplicate srcs (cdfs12220_broad.src): 
Input file with reference srcs (cdfs2313_broad.src): 
Transform file (out.xform): 

The output file, out.xform is as discussed above. The important parameters, the matrix elements look like

unix% dmlist out.xform"[cols a11,a12,a21,a22,t1,t2]" data,clean
#  a11                  a12                  a21                  a22                  t1                   t2
                  1.0                    0                    0                  1.0     0.32384497821605     0.31594220488116

This shows that to best match the locations between the two source lists, a shift 0.32 pixels in X and 0.32 pixels in Y direction is necessary. This will be accomplished with the wcs_update tool.



Apply transformation to Chandra Dataset

At this point in the thread the wcs_match output file out.xform has been created either by matching sources between a Chandra dataset and an external catalog or by matching sources between two Chandra datasets.

To update the Chandra dataset requires updating the aspect solution as well as the event file.

[WARNING]
Earlier versions of this thread

Earlier versions of this thread provided notes on correcting either the aspect solution or correcting the event file (and derived products). While those approaches worked for generating correct sky coordinates, the information needed to convert to chip coordinates and off-axis angles was lost.

The approach taken in this new version of the thread is the most robust way to update suite of Chandra data products for a given ObsID.

Update aspect solution values

First, the aspect solution file is updated. A new file, outfile, is created since every row in the input file will be modified.

unix% pset wcs_update infile=pcadf12220_000N001_asol1.fits
unix% pset wcs_update outfile=pcadf12220_000N001_corrected_asol1.fits
unix% pset wcs_update transformfile=out.xform
unix% pset wcs_update wcsfile=cdfs12220_broad_thresh.img
unix% wcs_update
Input coordinate transform file (out.xform): 
Either input asol file, or file with WCS to be updated (pcadf12220_000N001_asol1.fits): 
Output asol file (pcadf12220_000N001_corrected_asol1.fits): 
[IMPORTANT]
wcsfile

As was note earlier, make sure to use the same wcsfile as was used with wcs_match.

The dmdiff tool can be used to compare the difference in the aspect solution files.

unix% dmdiff pcadf12220_000N001_asol1.fits pcadf12220_000N001_corrected_asol1.fits  | less
...
RA_NOM    Values are not equal                     53.11785338186  53.117803347224  -5.00346e-05 (-9.42e-05%)
DEC_NOM   Values are not equal                   -27.802188221278  -27.802145042501 +4.31788e-05 (-0.000155%)
RA_PNT    Values are not equal                     53.11785338186  53.117803347224  -5.00346e-05 (-9.42e-05%)
DEC_PNT   Values are not equal                   -27.802188221278  -27.802145042501 +4.31788e-05 (-0.000155%)
...
----------------------------------------------------------------------
Compare Tables
----------------------------------------------------------------------
Compare Table Structure:
  Block name: ASPSOL
Compare Column Details:
Compare Virtual Column Details:
Compare Column Data:
Column:          Row:Cell:      Message:                                                        Value(s):                    Diff:
---------------- -------------- -------------------------------------- ---------------------------------- ------------------------
ra                      1       Values are not equal                    53.114386573516  53.1143365397926 -5.00337e-05 (-9.42e-05%)
dec                     1       Values are not equal                   -27.8031976810397  -27.8031545010137   +4.318e-05 (-0.000155%)
q_att                   1:1     Values are not equal                   0.798006784112701  0.798007118452388  +3.3434e-07 (+4.19e-05%)
q_att                   1:2     Values are not equal                   0.478068523024739  0.478068116791546 -4.06233e-07 (-8.5e-05%)
q_att                   1:3     Values are not equal                   -0.0773754694106698  -0.077375301446  +1.67965e-07 (-0.000217%)
q_att                   1:4     Values are not equal                   0.358676311646766  0.358676145475055 -1.66172e-07 (-4.63e-05%)
ra                      2       Values are not equal                   53.1143855993814  53.1143355656575 -5.00337e-05 (-9.42e-05%)
dec                     2       Values are not equal                   -27.8031998642721  -27.8031566842457   +4.318e-05 (-0.000155%)
q_att                   2:1     Values are not equal                   0.798006804947779  0.798007139287462  +3.3434e-07 (+4.19e-05%)
q_att                   2:2     Values are not equal                   0.478068514175583  0.478068107942399 -4.06233e-07 (-8.5e-05%)
q_att                   2:3     Values are not equal                   -0.0773755196886727  -0.0773753517239779 +1.67965e-07 (-0.000217%)
q_att                   2:4     Values are not equal                   0.358676266240039  0.358676100068302 -1.66172e-07 (-4.63e-05%)
...
[NOTE]
Changes in CIAO 4.8

In CIAO 4.8 both the RA_NOM & DEC_NOM and RA_PNT & DEC_PNT keywords in the aspect solution file are updated.

If there is more than 1 pcad_asol1.fits file, then each must be updated individually with wcs_update.


Update event file WCS

Next the World Coordinate System (WCS) parameters in the event file are updated. Since this only involves updating two header keyword values, the event file is modified in place. The outfile parameter should be set to "" (blank, no quotes). It is generally a good idea to make a copy of the event file and update the copy. Using dmcopy will copy the EVENTS extensions in the file and will provide a HISTORY record showing that file was copied.

unix% dmcopy acisf12220_repro_evt2.fits acisf12220_corrected_evt2.fits

unix% pset wcs_update infile=acisf12220_corrected_evt2.fits
unix% pset wcs_update outfile=
unix% wcs_update
Input coordinate transform file (out.xform): 
Either input asol file, or file with WCS to be updated (acisf12220_corrected_evt2.fits): 
Output asol file (): 
[IMPORTANT]
Avoid using dmcopy opt=all

Users with gratings data will need to be sure to also copy the REGION extension when copying the event file; however, due to a bug in dmcopy the use of the opt=all parameter is strongly discouraged. Instead uses should follow the suggested work around to use dmappend to copy the REGION block

unix% dmcopy evt.fits corrected_evt.fits
unix% dmappend "evt.fits[REGION][subspace -time]" corrected_evt.fits

Using dmdiff shows that there is a difference in the corrected file

unix% dmdiff acisf12220_repro_evt2.fits acisf12220_corrected_evt2.fits
...
RA_PNT   Values are not equal                     53.11785338186  53.117803347224  -5.00346e-05 (-9.42e-05%)
DEC_PNT  Values are not equal                   -27.802188221278  -27.802145042501 +4.31788e-05 (-0.000155%)
RA_NOM   Values are not equal                     53.11785338186  53.117803347224  -5.00346e-05 (-9.42e-05%)
DEC_NOM  Values are not equal                   -27.802188221278  -27.802145042501 +4.31788e-05 (-0.000155%)
...

----------------------------------------------------------------------
Compare Tables
----------------------------------------------------------------------
Compare Table Structure:
  Block name: EVENTS
Compare Column Details:
Compare Virtual Column Details:
# dmdiff (CIAO): WARNING: vector coordinate 'EQPOS' transform crval[0] differ.
# dmdiff (CIAO):   crval1=53.1179
# dmdiff (CIAO):   crval2=53.1178
# dmdiff (CIAO): WARNING: vector coordinate 'EQPOS' transform crval[1] differ.
# dmdiff (CIAO):   crval1=-27.8022
# dmdiff (CIAO):   crval2=-27.8021
Compare Column Data:

however the precision is insufficient to see the actual numeric difference. The raw keywords can be inspected, specifically looking for the TCRVLnn keywords

unix% dmlist acisf12220_repro_evt2.fits header,clean,raw | grep CRVL
TCRVL5       = 0                    /                     
TCRVL6       = 0                    /                     
TCRVL9       = 0                    /                     
TCRVL10      = 0                    /                     
TCRVL11      =      53.11785338     /                     
TCRVL12      =     -27.80218822     /                     
unix% dmlist acisf12220_corrected_evt2.fits header,clean,raw | grep CRVL
TCRVL5       = 0                    /                     
TCRVL6       = 0                    /                     
TCRVL9       = 0                    /                     
TCRVL10      = 0                    /                     
TCRVL11      =      53.11780335     /                     
TCRVL12      =     -27.80214504     /                     
[NOTE]
Changes in CIAO 4.8

In CIAO 4.8 both the RA_NOM & DEC_NOM and RA_PNT & DEC_PNT keywords in the event file are updated.

To allow downstream analysis tasks to work correctly, users should update the ASOLFILE keyword in the updated event file with the updated aspect solution file name

unix% dmhedit acisf12220_corrected_evt2.fits file= op=add key=ASOLFILE value=pcadf12220_000N001_corrected_asol1.fits

This completes all the edits necessary to correct the fine astrometry for this Chandra dataset. At this point there is no need to run chandra_repro or any of the event processing tools (acis_process_events nor hrc_process_events). However, any data products that were derived from the event file before the correction need to be recreated. This will include the field of view files, images, and exposure maps.

Individual files may be updated using the same recipe as for the event file. For example:

unix% dmcopy cdfs2313_broad_flux.img cdfs2313_broad_flux.corrected.img op=all
unix% wcs_update cdfs2313_broad_flux.corrected.img outfile= trans=out.xform wcsfile=cdfs2313_broad_flux.img

Remember: the wcsfile must be the original, uncorrected image since that is what was used to generate the transformation coefficients. However, given the number of files it is likely easier to simply recreate them from the correct aspect solution and event file.


Alternative: Using offsets

It may be the case that users have a source in the field with a known source location. It is also possible to adjust the astrometry of the Chandra observation to match that location.

In this example the source shown in Figure 4 is assumed to be a known source with coordinates

03:32:38.0 -27:52:13.0

These coordinates are just an example for the purposes of running this thread.

Figure 4: Source with known coordinates

[Thumbnail image: Example source]

[Version: full-size]

[Print media version: Example source]

Figure 4: Source with known coordinates

The broad-band image for OBSID 12220 centered on a source whose coordinates are assumed to be known.

The offset between a source's measured location and the known location can also be input to wcs_update. First the source location is measured. It could have been identifed from the wavdetect source list used above; but for this example dmstat will be used to compute the source centroid

unix% dmstat cdfs12220_broad_thresh.img"[sky=circle(3840,3599,7)]" cen+ sig- med-
EVENTS_IMAGE(x, y)
    min:	0 	      @:	( 3798 3605 )
    max:	6 	      @:	( 3800 3613 )
cntrd[log] :	( 8.75 8.1166666667  )
cntrd[phys]:	( 3840.75 3599.1166667 )
   good:	149 
   null:	76 

The values of interest are the cntrd[phys] (centroid in physical coordinates) values. These can be used with dmcoords to convert the physical coordinates to RA/Dec

unix% punlearn dmcoords
unix% dmcoords cdfs12220_broad_thresh.img op=sky x=3840.75 y=3599.11 verb=0 celfmt=deg 
unix% pget dmcoords ra dec
53.15739196468656
-27.87015918504676

and then the prop_precess tool to is used to convert between degrees and sexagesimal format

unix% prop_precess f j/deg t j/hms eval 53.15739196468656 -27.87015918504676
Precess [Conversion mode]
Enter "q" to return to setup mode
--------------------------------------------------------------------------------
RA,Dec J2000.0                   53.157392   -27.870159
RA,Dec J2000.0                  03 32 37.77 -27 52 12.57
--------------------------------------------------------------------------------

This shows that observed source location is at 03:32:37.77 -27:52:12.57. dmcoords could have been run with with celfmt=hms to get the coordinates in sexagesimal format directly.

Next, the known source location, 03:32:38.0 -27:52:13.0, needs to be converted to physical coordinates . This is done directly with dmcoords using the celfmt=hms option:

unix% punlearn dmcoords
unix% dmcoords cdfs12220_broad_thresh.img op=cel celfmt=hms ra=03:32:38.0 dec=-27:52:13.0 verb=0
unix% pget dmcoords x y
3834.661168246338
3598.240259287044

The coordinates are also converted to decimal degrees using prop_precess

unix% prop_precess from j/hms to j/deg p0 eval 03 32 38.0 -27 52 13.0
53.158333   -27.870278

The prop_precess p0 option has been used which limits the amount of screen output. Note: prop_precess requires spaces, not colons (":") between the coordinate values.

The necessary inputs needed to compute the astrometric correction have been computed. As discussed above, the wcs_update tool uses physical pixel offsets defined relative to the input tangent plane (the wscfile). Specifically the offset from the reference coordinates to the measured coordinates needs to be computed. In this example

unix% echo 3834.661168246338 - 3840.75  | bc -l
-6.088831753662
unix% echo 3598.240259287044 - 3599.1166667 | bc -l
-.876407412956

These offsets can then be applied with wcs_update directly using the deltax and deltay parameters:

unix% dmcopy cdfs12220_broad_thresh.img cdfs12220_broad_thresh.img.manual_offsets clob+
unix% wcs_update cdfs12220_broad_thresh.img.manual_offsets none trans= wcsfile=")infile" deltax=-6.088831753662 deltay=-.876407412956

To check that the update was applied the dmstat and dmcoords commands are repeated on the updated output file:

unix% dmstat cdfs12220_broad_thresh.img"[sky=circle(3840,3599,7)]" cen+ sig- med-
EVENTS_IMAGE(x, y)
    min:	0 	      @:	( 3840 3592 )
    max:	6 	      @:	( 3841 3598 )
cntrd[log] :	( 8.75 8.1166666667 )
cntrd[phys]:	( 3840.75 3599.1166667 )
   good:	149 
   null:	76 

As expect the exact same centroid in physical coordinates is obtained as in the original image. The expected difference is in the world coordinates, RA/Dec output from dmcoords

unix% punlearn dmcoords
unix% dmcoords cdfs12220_broad_thresh.img.manual op=sky x=3840.75 y=3599.1166667 verb=0 celfmt=hms
unix% pget dmcoords ra dec
03:32:37.999
-27:52:12.99

which now matches the known location.

Why compute offsets by hand?

Rather than compute the offsets by hand, it is also easy to create a pair of ASCII files that can be input to wcs_match using the RA and Dec values in decimal degrees.

First, the reference coordinates

unix% cat << EOM > usr.dat
#ra dec
53.158333   -27.870278
EOM

and then the measured centroid coordinates

unix% cat << EOM > obs.dat
#ra dec
53.15739196468656 -27.87015918504676
EOM

and then run wcs_match

unix% wcs_match in=obs.dat ref=usr.dat out=manual.xform wcs=cdfs12220_broad_thresh.img  clob+ method=trans verb=1
 input (dup) src file : obs.dat
 input ref src file   : usr.dat
 input wcsfile        : cdfs12220_broad_thresh.img
 debug level          : 1

1 common sources found between: 
usr.dat
obs.dat
After deleting poor matches, 1 sources remain

Source Residuals
----------------
 Match Ref# Dup#    Ref RA      Ref Dec.    Prior Resid           Transfm Resid         Resid  Incl
 Index              (deg.)      (deg.)      RSS (x,y)             RSS (x,y)             Ratio
                                            (arcsec)              (arcsec)
   0    0     0     53.15833   -27.87028    3.03 ( 2.99, 0.43)    0.00 (-0.00,-0.00)    0.00    Y

Source Residuals, before/after transform (arcsec), and percentage improvement:

   Average Residuals:         3.025176   0.000000  100.00%
   Maximum Residuals:         3.025176   0.000000  100.00%
   RMS Residuals:             2.139122   0.000000  100.00%

Source Residual Ratios, before/after transform, and percentage improvement:

   Average Residual Ratios:   3.074365   0.000000  100.00%
   Maximum Residual Ratios:   3.074365   0.000000  100.00%
   RMS Ratios:                2.173905   0.000000  100.00%

# wcs_match (CIAO): WARNING: Dup src x_err or ra_err cols not found.  Assuming x_err = 1.
# wcs_match (CIAO): WARNING: Dup src y_err or dec_err cols not found.  Assuming y_err = 1.
# wcs_match (CIAO): WARNING: Ref src x_err or ra_err cols not found.  Assuming x_err = 1.
# wcs_match (CIAO): WARNING: Ref src y_err or dec_err cols not found.  Assuming y_err = 1.

These ERROR and WARNING messages can be ignored here since the method=trans option is used which only requires a single matching source pair. The output file transform file looks like

unix% dmlist manual.xform data,clean 
#  a11                  a12                  a21                  a22                  t1                   t2                   ra_ref               dec_ref              roll_ref             xpix_ref             ypix_ref             x_scale              y_scale
                  1.0                    0                    0                  1.0        -6.0866750980    -0.87136601986367        53.1178533819       -27.8021882213                    0              4096.50              4096.50     -0.0001366666667      0.0001366666667

The offset columns: t1 and t2 show the same values that were computed by hand above. This file can then be used in wcs_update

unix% dmcopy cdfs12220_broad_thresh.img cdfs12220_broad_thresh.img.manual clob+
unix% wcs_update cdfs12220_broad_thresh.img.manual none trans=manual.xform wcsfile=")infile" 

and the same coordinate verifications performed:

unix% dmstat cdfs12220_broad_thresh.img"[sky=circle(3798,3612,7)]" cen+ sig- med-
EVENTS_IMAGE(x, y)
    min:	0 	      @:	( 3840 3592 )
    max:	6 	      @:	( 3841 3598 )
cntrd[log] :	( 8.75 8.1166666667 )
cntrd[phys]:	( 3840.75 3599.1166667 )
   good:	149 
   null:	76 

unix% punlearn dmcoords
unix% dmcoords cdfs12220_broad_thresh.img.manual op=sky x=3840.75 y=3599.1166667 verb=0 celfmt=hms
unix% pget dmcoords ra dec
03:32:37.999
-27:52:12.99

Which again confirms the correct offsets have been applied.


Summary

This thread has shown how to apply fine astrometic shifts to a Chandra dataset.

This involves generating two source lists: one for the current dataset and a reference list that could be either a Chandra dataset or some external catalog. While wavdetect was use here, any source detection algorithm may be used.

Then a transformation matrix is derived by cross matching the sources using wcs_match.

This transformation matrix is then applied to the aspect solution file and the event file using wcs_update.


History

24 Apr 2006 new for CIAO 3.3: original version
01 Dec 2006 updated for CIAO 3.4: CIAO version in error messages
23 Jan 2008 updated for CIAO 4.0: kernel parameter removed from wavdetect; filename and screen output updated for reprocessed data (version N003 event file for 441)
23 Jun 2008 updated image display to place figures inline with text
13 Feb 2009 updated for CIAO 4.1: minor change to screen output due to CIAO 4.1 change in reproject_aspect
22 May 2009 added a note on combining the data.
22 Jul 2009 clarified that the source list does not have to be made by wavdetect; a list from another data analysis tool or catalog may be used
08 Feb 2010 reviewed for CIAO 4.2: no changes
13 Jan 2011 reviewed for CIAO 4.3: no changes
11 Jan 2012 reviewed for CIAO 4.4: minor change to output from wavdetect (34 common sources found by reproject_aspect instead of 33)
03 Dec 2012 Review for CIAO 4.5; no changes
02 May 2013 Updated with example showing how to use new asol file with chandra_repro; clarified allowed format (FITS or ASCII); added example of filtering source lists;
04 Apr 2014 Added using mkpsfmap for the psffile parameter used by wavdetect.
28 Jul 2014 This thread has been completely rewritten. This version corrects a problem with the reprojected data not properly updating the RA_PNT and DEC_PNT keywords which can lead to incorrect coordinate transformations. The thread now uses wcs_update and wcs_match separately since the former needs to be run multiple times. It also shows how to do a match with an optical catalog.
15 Dec 2014 Reviewed for CIAO 4.7. The change in wcs_update to modify the *_NOM keywords simplifies the end of this thread.
23 Nov 2015 Reviewed for CIAO 4.8. Removed extra steps to update the _NOM and _PNT keywords, wcs_update does that automatically now.
14 Jul 2017 Added a new section showing how to use wcs_update offset parameters directly.
13 May 2019 Added note about new fluximage psfecf parameter in CIAO 4.11.
25 Jan 2021 Added link to absolute astrometric accuracy page.
22 Feb 2022 Review for CIAO 4.14. Updated for Repro5/CALDB 4.9.6