[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.
>
> 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

I used astrometry.net to get the WCS headers, and for that I had to crop 
the original image such that only the starfield is visible (the original 
image contained some 'earthly' structures). But after that I actually 
need coordinates of the whole original image, that is, pixels which were 
outside the image used for astrometry solving. So it sounds like your 
latter change should be useful for that case.

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