As soon as I become more and more observational astronomer, I experience complicated problems such as correcting astrometry of objects in the XMM-Newton field of view. The source positions provided by the telescope are already very good quality and presumably do not differ from actual positions by more than 3.6 arcsec. Nevertheless, there are many cases when the further improvement is desirable. If the Chandra observatory is used the standard package CIAO already has tools to solve this problem, they are reproject_aspect and wcs_match. The standard software of the XMM-Newton (SAS) has the only a tool to recalibrate astrometry based on the optical image obtained by the telescope itself. Nothing is written in manuals about astrometry correction using a list of sources from an external catalog such as the USNO or 2MASS.
It appears that this problem can be solved by means of python and amazing astropy package. The following links were the most useful to understand the principles of WCS (World coordinate system) in FITS (flexible image transport system): old blog post with using pywcs, stack exchange post applied to NuStar observations and example of how to use astropy.wcs.
The idea behind the script is quite simple. The astropy library has class wcs which can be initialized directly using the fits header. The wcs object represents a transformation of pixel number on CCD array to right ascension and declination including optical distortions and many more tiny instrumental effects. The information in wcs might be slightly wrong. If you have an astrometric catalog of bright sources, you can improve this information slightly changing two wcs parameters: CRVAL and PC value. First one – is the coordinate reference value and is a vector. The second value is a rotation and is represented by a rotational matrix. After we change both these values we compute new sky coordinate using function all_pix2world. These new coordinates are compared with ones obtained from a catalog and residuals are optimized using e.g. least_squares from scipy.optimize.
Here is an example of the code:
from scipy.optimize import least_squares
from math import *
import numpy as np
import csv
import pywcs
import pyfits
from astropy.wcs import WCS
import copy
import sys
w = None
## Function which describes a transformation from pixel coordinate to the right ascension and declination
def tranf (param, x_vect):
phi = param[0] ## Rotation angle
A1 = param[1] ## Shift
A2 = param[2]
original_crval = w.wcs.crval.copy() ## Keep original value
original_pc = w.wcs.pc.copy() ## here
w.wcs.crval = original_crval.copy() + [A1, A2]
to_matr = [[cos(phi), sin(phi)], [-sin(phi), cos(phi)]] ## The original matrix was identical here, but in principle it should be a matrix multiplication.
w.wcs.pc = to_matr
ra_pred, dec_pred = w.all_pix2world(x_vect[0], x_vect[1], 0) ## Compute the sky coordinate. The last value means that we start counting pixels from 0.
w.wcs.crval = original_crval.copy() ## Restore the original values
w.wcs.pc = original_pc.copy()
return [ra_pred, dec_pred]
## Here we compute residuals between coordinates in catalogue (y) and coordinates derived from the pixel number (x) using parameters of transfromation param
def resid (param, x, y):
res = []
for i in range (0, len(x)):
ra_pred, dec_pred = tranf (param, x[i])
d_ra = (ra_pred - y[i][0]) * 3600 ## The residuals are in arcsec
d_dec = (dec_pred - y[i][1]) * 3600
#print 'RA, DEC (predicted, actual): ', ra_pred, dec_pred, ' \t ', y[i][0], y[i][1]
resid = sqrt(pow(d_ra,2) + pow(d_dec,2.0))
res.append(resid)
return res
x=[]
y=[]
ra=[]
dec=[]
counter = 0
with open('catalog.csv', 'rb') as f: ## Cross-match between sources detected in my field by the XMM-Newton and USNO B1
reader = csv.reader(f)
for row in reader:
if counter == 0:
for i in range (0, len(row)): ## Check that is the header of the table
print i, row[i]
counter = counter + 1
else:
if row[13] != '' and row[14]!= '':
x.append (float(row[2])) ## coordinates at the CCD
y.append (float(row[3]))
ra.append (float(row[8])) ## coordinates at the sky
dec.append(float(row[9]))
N = len(x)
w = WCS('epn-s.fits') ## Read WCS information from fits file
x0 = [0, 0, 0] ## Initial guess
before_vect = np.zeros((N,2))
after_vect = np.zeros((N,2))
for i in range (0, N): ## Initialise the x,y coordinates at CCD and alpha, delta at the sky
before_vect[i][0] = x[i]
before_vect[i][1] = y[i]
after_vect[i][0] = ra[i]
after_vect[i][1] = dec[i]
res_lsq = least_squares(resid, x0, args=(before_vect, after_vect), max_nfev=30000, method='lm', xtol=1e-15, ftol=1e-15) ## Optimisation
print '--------------'
print res_lsq
final = res_lsq.x
val= resid (final, before_vect, after_vect)
print 'Std: ', np.std(val) ## Check how good our fit is