[AstroPy] Running wcs_pix2world for all image pixels

Fri Oct 18 04:54:16 EDT 2013

On 17 October 2013 16:28, Maik Riechert <maik.riechert at arcor.de> wrote:

I didn't use tranGrid here, just tran! So no tolerance involved.
>

Doh! Sorry - didn't spot that.

>
> NAXIS1  =                 1468 / image
> width
> NAXIS2  =                  916 / image
> height
> CTYPE1  = 'RA---TAN-SIP' / TAN (gnomic) projection + SIP distortions
> CTYPE2  = 'DEC--TAN-SIP' / TAN (gnomic) projection + SIP distortions
> EQUINOX =               2000.0 / Equatorial coordinates definition (yr)
> LONPOLE =                180.0 / no
> comment
> LATPOLE =                  0.0 / no
> comment
> CRVAL1  =         302.88771538 / RA  of reference
> point
> CRVAL2  =        41.8941029213 / DEC of reference
> point
> CRPIX1  =                  735 / X reference
> pixel
> CRPIX2  =                  459 / Y reference
> pixel
> CUNIT1  = 'deg     ' / X pixel scale
> units
> CUNIT2  = 'deg     ' / Y pixel scale
> units
> CD1_1   =     -0.0157184611387 / Transformation
> matrix
> CD1_2   =     0.00022032143241 / no
> comment
> CD2_1   =   -0.000284856284919 / no
> comment
> CD2_2   =     -0.0159943489114 / no
> comment
> IMAGEW  =                 1468 / Image width,  in
> pixels.
> IMAGEH  =                  916 / Image height, in
> pixels.
> A_ORDER =                    2 / Polynomial order, axis
> 1
> A_0_2   =   -5.10167355829E-05 / no
> comment
> A_1_1   =   -5.54267287337E-05 / no
> comment
> A_2_0   =    3.68095680847E-05 / no
> comment
> B_ORDER =                    2 / Polynomial order, axis
> 2
> B_0_2   =    2.49853669312E-05 / no
> comment
> B_1_1   =    6.07756179822E-05 / no
> comment
> B_2_0   =   -1.70350503136E-05 / no
> comment
> AP_ORDER=                    2 / Inv polynomial order, axis
> 1
> AP_0_1  =     -0.0014471978634 / no
> comment
> AP_0_2  =     5.0935282476E-05 / no
> comment
> AP_1_0  =    0.000641040653934 / no
> comment
> AP_1_1  =    5.54267224252E-05 / no
> comment
> AP_2_0  =     -3.674588665E-05 / no
> comment
> BP_ORDER=                    2 / Inv polynomial order, axis
> 2
> BP_0_1  =      0.0010860421788 / no
> comment
> BP_0_2  =   -2.49999638422E-05 / no
> comment
> BP_1_0  =   -0.000569392740582 / no
> comment
> BP_1_1  =   -6.07473585161E-05 / no comment
> BP_2_0  =     1.7000827442E-05 / no comment
>

OK - so I suspect the problems may actually be in astropy rather than AST.
 Maybe the maintainers of astropy.wcs could comment.

The above header includes a polynomial distortion correction using the
Spitzer SIP conventions. The "A_" and "B_" headers give the coefficients of
the forward ( (x,y) to (ra,dec) ) polynomial, and the "AP_" and "BP_"
headers give the coefficients of the inverse ( (ra,dec) to (x,y) )
polynomial.

Using the header as supplied to transform the following pixel positions:

x = -196, 735, 4059
y = -504, 459, 2328

I get, using AST (converted to degrees):

ra = 328.090902, 302.887715, 265.839543
dec = 53.570723, 41.894103, 8.257963

whereas astropy gives:

ra = 326.992376, 302.866900, 263.545072
dec = 54.749009, 41.877821, 10.589047

For an independent check I also used SAO WCSTools which gave:

ra =  328.09090, 302.88772, 265.83954
dec = 53.57072, 41.89410, 8.25796

So AST and WCSTools agree, but astropy is different. This is not conclusive
in itself, but I then modified the header by removing all the SIP keywords
and removing "-SIP" from the CTYPE1/2 headers. With these changes the
header represents an undistorted TAN projection. I then transformed the
same pixel positions again using the changed header:

AST gave:

ra = 327.022569, 302.887715, 263.548641
dec = 54.758313, 41.894103, 10.599558

astropy gave:

ra = 326.992376, 302.866900, 263.545072
dec = 54.749009, 41.877822, 10.589047

and WCSTools gave:

ra = 327.02257, 302.88772, 263.54864
dec = 54.75831, 41.89410, 10.59956

So again AST and WCSTools agree, with astropy as the odd one out. But more
importantly, the astropy results are exactly the same as when using the
header that includes the SIP distortion. This seems to indicate that in
fact astropy is ignoring the SIP distortion. You can see that the astropy
values are much closer to the AST/WCSTools values when no distortion is
included.

Further, the second pixel position is actually the reference point of the
header (CRPIX1/2) and so the RA and Dec at this point should be exactly
equal to the values of CRVAL1/2, which is what AST and WCSTools report. But
astropy is off by about a pixel (58 arc-sec).

So to summarise, astropy seems to ignore the SIP distortion, and also
assigns a slightly wrong position to the reference point.  The astropy code
I am using is:

from astropy import wcs
from astropy.io import fits
hdulist = fits.open("maik.fit")
mywcs = wcs.WCS(hdulist[0].header)
x = [ -196.,   735.,  4059.]
y = [ -504.,   459.,  2328.]
ra, dec = mywcs.wcs_pix2world(x, y, 0)

Now regarding the inverse transformation back from (ra,dec) to (x,y). If in
fact astropy is not including any polynomial distortion, then it is
not surprising that it manages to recreate the original (x,y) values
exactly. AST also does that if you exclude the distortion terms. If the
distortion is included, AST gives:

x = -205.93, 735, 4179.93
y = -499.37, 459, 2246.86

WCSTools gives exactly the same values. This poor performance of AST and
WCSTools when including the distortion indicates that the coefficients of
the inverse polynomial included in the header (the "AP_" and "BP_" headers)
are inaccurate.  But AST has a facility to handle SIP headers that do not
include an inverse polynomial. In such cases, AST creates an inverse itself
by numerical inversion of the supplied forward polynomial. So if I remove
the "AP_" and "BP_" keywords from the header, and again transform the
(ra,dec) values back into (x,y) AST now gives:

x = -197.19, 735, 4061.94
y = -502.92, 459, 2324.94

which are much better but still not good. AST includes two alternative
schemes for implementing a missing SIP inverse:

1) it samples the forward transformation at a regular grid of pixel
positions and fits a polynomial surface to the results. This fit is
performed once, and if the residuals of the fit are sufficiently small, it
is then used to transform all subsequent points.

2) if the residuals of the above fit are too large, then each individual
(ra,dec) position is transformed using an iterative Newton-Raphson method.
This is much slower over-all than 1) but more accurate.

There is a problem in AST in that it is using option 1) above for your
header, when in fact it needs to use 2) to get sufficient accuracy.  The
recovered (x,y) values quoted above were created using option 1). With
option 2), the exact original (x,y) values are recovered to within 0.01
pixel.

I shall investigate this issue in AST. Thanks for the input.

David
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/astropy/attachments/20131018/3d90291b/attachment.html>