{"id":43,"date":"2022-02-18T15:28:14","date_gmt":"2022-02-18T15:28:14","guid":{"rendered":"http:\/\/www.pulsars.info\/wordpress\/?p=43"},"modified":"2022-02-18T15:28:37","modified_gmt":"2022-02-18T15:28:37","slug":"correcting-astrometry","status":"publish","type":"post","link":"http:\/\/www.pulsars.info\/wordpress\/uncategorized\/correcting-astrometry\/","title":{"rendered":"Correcting astrometry"},"content":{"rendered":"\n
\"\"<\/figure><\/div>\n\n\n\n

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.<\/p>\n\n\n\n

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<\/a>, stack exchange post applied to NuStar observations<\/a> and example of how to use astropy.wcs<\/a>.<\/p>\n\n\n\n

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.<\/p>\n\n\n\n

Here is an example of the code:<\/p>\n\n\n\n

from<\/strong> scipy.optimize import<\/strong> least_squares\nfrom<\/strong> math import<\/strong> *\nimport<\/strong> numpy as<\/strong> np\nimport<\/strong> csv\nimport<\/strong> pywcs\nimport<\/strong> pyfits\nfrom<\/strong> astropy.wcs import<\/strong> WCS\nimport<\/strong> copy\nimport<\/strong> sys\n\nw = None\n\n## Function which describes a transformation from pixel coordinate to the right ascension and declination<\/em>\n\ndef<\/strong> tranf (param, x_vect):\n\n    phi = param[0] ## Rotation angle<\/em>\n    A1  = param[1] ## Shift <\/em>\n    A2  = param[2]\n\n    original_crval = w.wcs.crval.copy() ## Keep original value<\/em>\n    original_pc    = w.wcs.pc.copy()    ## here<\/em>\n\n    w.wcs.crval = original_crval.copy() + [A1, A2]\n    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.<\/em>\n    w.wcs.pc = to_matr\n\n    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.<\/em>\n\n    w.wcs.crval = original_crval.copy() ## Restore the original values <\/em>\n    w.wcs.pc    = original_pc.copy()\n\n    return<\/strong> [ra_pred, dec_pred]\n\n## Here we compute residuals between coordinates in catalogue (y) and coordinates derived from the pixel number (x) using parameters of transfromation param<\/em>\ndef<\/strong> resid (param, x, y):\n\n    res = []\n    for<\/strong> i in<\/strong> range (0, len(x)):\n\n        ra_pred, dec_pred = tranf (param, x[i]) \n        d_ra  = (ra_pred  - y[i][0]) * 3600  ## The residuals are in arcsec<\/em>\n        d_dec = (dec_pred - y[i][1]) * 3600\n\n        #print 'RA, DEC (predicted, actual): ', ra_pred, dec_pred, ' \\t ', y[i][0], y[i][1]<\/em>\n\n        resid = sqrt(pow(d_ra,2) + pow(d_dec,2.0))\n        res.append(resid)\n\n    return<\/strong> res\n\n\nx=[]\ny=[]\nra=[]\ndec=[]\n\ncounter = 0\n\n\nwith<\/strong> open('catalog.csv', 'rb') as<\/strong> f:  ## Cross-match between sources detected in my field by the XMM-Newton and USNO B1<\/em>\n    reader = csv.reader(f)\n    for<\/strong> row in<\/strong> reader:\n\n    if<\/strong> counter == 0:\n        for<\/strong> i in<\/strong> range (0, len(row)): ## Check that is the header of the table<\/em>\n            print<\/strong> i, row[i]\n        counter = counter + 1\n\n    else<\/strong>:\n\n        if<\/strong> row[13] != '' and<\/strong> row[14]!= '':\n            x.append  (float(row[2])) ## coordinates at the CCD<\/em>\n                        y.append  (float(row[3]))\n            ra.append (float(row[8])) ## coordinates at the sky<\/em>\n            dec.append(float(row[9]))\n\n\nN = len(x)\n\nw = WCS('epn-s.fits') ## Read WCS information from fits file<\/em>\n\nx0 = [0,  0, 0] ## Initial guess<\/em>\n\nbefore_vect = np.zeros((N,2))\nafter_vect  = np.zeros((N,2))\n\nfor<\/strong> i in<\/strong> range (0, N):  ## Initialise the x,y coordinates at CCD and alpha, delta at the sky<\/em>\n    before_vect[i][0] = x[i]\n    before_vect[i][1] = y[i]\n    after_vect[i][0] = ra[i]\n    after_vect[i][1] = dec[i]\n\nres_lsq = least_squares(resid, x0, args=(before_vect, after_vect), max_nfev=30000, method='lm', xtol=1e-15, ftol=1e-15) ## Optimisation<\/em>\n\nprint<\/strong> '--------------'\nprint<\/strong> res_lsq\n\nfinal = res_lsq.x\n\nval= resid (final, before_vect, after_vect)\n\nprint<\/strong> 'Std: ', np.std(val)  ## Check how good our fit is <\/em>\n<\/pre>\n","protected":false},"excerpt":{"rendered":"
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…<\/p>\n
Continue ReadingCorrecting astrometry<\/span><\/a><\/div>\n","protected":false},"author":1,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[1],"tags":[],"_links":{"self":[{"href":"http:\/\/www.pulsars.info\/wordpress\/wp-json\/wp\/v2\/posts\/43"}],"collection":[{"href":"http:\/\/www.pulsars.info\/wordpress\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"http:\/\/www.pulsars.info\/wordpress\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"http:\/\/www.pulsars.info\/wordpress\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"http:\/\/www.pulsars.info\/wordpress\/wp-json\/wp\/v2\/comments?post=43"}],"version-history":[{"count":1,"href":"http:\/\/www.pulsars.info\/wordpress\/wp-json\/wp\/v2\/posts\/43\/revisions"}],"predecessor-version":[{"id":44,"href":"http:\/\/www.pulsars.info\/wordpress\/wp-json\/wp\/v2\/posts\/43\/revisions\/44"}],"wp:attachment":[{"href":"http:\/\/www.pulsars.info\/wordpress\/wp-json\/wp\/v2\/media?parent=43"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"http:\/\/www.pulsars.info\/wordpress\/wp-json\/wp\/v2\/categories?post=43"},{"taxonomy":"post_tag","embeddable":true,"href":"http:\/\/www.pulsars.info\/wordpress\/wp-json\/wp\/v2\/tags?post=43"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}