[AstroPy] Running wcs_pix2world for all image pixels

[Maik Riechert]

On 17/10/13 08:54, David Berry wrote:
> AST attempts to speed up this sort of thing by using a recursive 
> scheme as follows:
> ...
> This means that for typical mappings that are locally smooth, you are 
> usually using a very simple and fast linear transformation for each 
> point rather than the full might of the complete pixel->wcs mapping. 
> Clearly, you need to set some limit to the depth of recursion - we use 
> the criterion that the number of pixels in the pane must be at least 4 
> times the number of points needed to do the linear fit.
>
> The user specifies the accuracy with which the points are to be found, 
> and this determines the maximum acceptable residuals for the fit. If 
> the residuals are higher, the pane is subdivided.
>
> See for instance 
> http://starlink.jach.hawaii.edu/devdocs/sun211.htx/node430.html

I tried to use trangrid (with starlink-pyast) but I must be doing 
something wrong:

width = hdulist[0].header['NAXIS1']
height = hdulist[0].header['NAXIS2']
(wcs,_) = starlink.Atl.readfitswcs(hdulist[0])
mapping = wcs.getmapping( starlink.Ast.BASE, starlink.Ast.CURRENT )
print mapping.trangrid([0,0],[width-1,height-1])

[[-0.69877269 -0.69914315 -0.69951373 ..., -1.22724324 -1.22756251
   -1.22788175]
  [ 0.83534859  0.83540118  0.83545372 ...,  0.5841799   0.58413968
    0.58409942]]

I am expecting coordinates ranging from 320..285 for RA and 48..33 for Dec.

For comparison, my astropy code looks like that:

wcs = WCS(hdulist[0].header)
x, y = np.meshgrid(np.arange(width), np.arange(height))
ra, dec = wcs.wcs_pix2world(x, y, 0)

Cheers
Maik

> David
>
>     Cheers,
>     Tom
>
>     >
>     > Maik
>     >
>     >
>     > On 16/10/13 15:15, Thomas Robitaille wrote:
>     >>
>     >> Hi Maik,
>     >>
>     >> Someone asked a very similar question just yesterday on GitHub:
>     >>
>     >> https://github.com/astropy/astropy/issues/1587
>     >>
>     >> In short, the solution in your case is:
>     >>
>     >> NAXIS1 = 4256
>     >> NAXIS2 = 2832
>     >> x = np.arange(NAXIS1)
>     >> y = np.arange(NAXIS2)
>     >> X, Y = np.meshgrid(x, y)
>     >> ra, dec = wcs.wcs_pix2world(X, Y, 0)
>     >>
>     >> Any luck?
>     >>
>     >> Tom
>     >>
>     >>
>     >>
>     >>
>     >> On 16 October 2013 14:48, Maik Riechert <maik.riechert at arcor.de
>     <mailto:maik.riechert at arcor.de>> wrote:
>     >>>
>     >>> Hi,
>     >>>
>     >>> I'm quite new to python and astropy and got stuck at the
>     following which
>     >>> probably is extremely easy to do.
>     >>>
>     >>> I created a WCS object out of a FITS header and wanted to get
>     the RA/Dec
>     >>> coordinates of each image pixel. I tried:
>     >>>
>     >>> w.wcs_pix2world(np.ndindex(4256, 2832), 0)
>     >>>
>     >>> where 4256x2832 are the image dimensions. This returns
>     'ValueError:
>     >>> object of too small depth for desired array'. I was trying to
>     avoid
>     >>> allocating an array with all pixel coordinates and thought an
>     iterator
>     >>> would work too.
>     >>>
>     >>> What would you recommend for my case (translating every pixel)?
>     >>>
>     >>> Cheers and thanks,
>     >>> Maik
>     >>>
>     >>>
>     >>> _______________________________________________
>     >>> AstroPy mailing list
>     >>> AstroPy at scipy.org <mailto:AstroPy at scipy.org>
>     >>> http://mail.scipy.org/mailman/listinfo/astropy
>     >
>     >
>     _______________________________________________
>     AstroPy mailing list
>     AstroPy at scipy.org <mailto:AstroPy at scipy.org>
>     http://mail.scipy.org/mailman/listinfo/astropy
>
>

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