[AstroPy] Running wcs_pix2world for all image pixels

[David Berry]

On 18 October 2013 09:54, David Berry <d.berry at jach.hawaii.edu> wrote:

>
> Now regarding the inverse transformation back from (ra,dec) to (x,y).
> [...] 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.
>

In case you are interested, I've modified AST (v7.3.3) and pyast (v2.5) to
fix this.

Since I have come across several SIP headers that contain badly inaccurate
values for the inverse polynomial, I've followed astropy's lead and changed
AST so that it *always* calculates it's own inverse, ignoring any inverse
in the header.

The errors in the reconstructed (x,y) values were caused by the points
being way outside the bounds of the image. When AST calculates a fit to the
inverse, it previously used a grid that had the same bounds as the image,
so points outside this area were under-constrained, leading to poor fits
and bad errors. I've changed it so that the fit is now done on a much
larger area, and to a much  higher accuracy.  In the particular case of
your header, the fit fails to achieve the required accuracy and so the
transformation reverts to the Newton Raphson approach, which is slower but
more accurate.

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