[AstroPy] Running wcs_pix2world for all image pixels

[Maik Riechert]

On 18/10/13 14:48, David Berry wrote:
>
> On 18 October 2013 09:54, David Berry <d.berry at jach.hawaii.edu 
> <mailto: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.

Just a final remark. Having changed the origin to 1 for astropy I could 
produce consistent RA,Dec results from astropy (with SIP), starlink, and 
esutil. For starlink, the reverse direction is consistent now too, with 
the new version. A quick measurement yielded 3.6s for AstroPy, 3.8s for 
starlink (not using approximations), and 6.3s for esutil.

Thanks again for your support, it helped a lot!
Maik
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/astropy/attachments/20131021/2980e371/attachment.html>