[AstroPy] Running wcs_pix2world for all image pixels

[David Berry]

On 16 October 2013 22:02, Thomas Robitaille <thomas.robitaille at gmail.com>wrote:

> On 16 October 2013 15:58, Maik Riechert <maik.riechert at arcor.de> wrote:
> > Hi Tom,
> >
> > That's exactly what I needed. Thanks!
> >
> > It takes around 3s on my machine for a 4200x2800 image. Is that to be
> > expected or are there some speedup tricks?
>
> To be fair, you *are* converting 12 million positions through a WCS
> transformation which is a non-trivial transformation so I'm not sure
> if there are any speedups if you need all the values. However, the
> question is whether you really need the world coordinates for *all*
> pixels. What are you using the values for?
>
>

AST attempts to speed up this sort of thing by using a recursive scheme as
follows:



Initial pane = whole array
wcs_values = TransformPane( Mapping, pane )


where the TransformPane function does the following:



- Do a linear fit between pixel and world coords at a small number of
positions spread across the pane.
- If the residuals of the fit are small enough:
   - use the fit to transform all positions in the pane and return them as
the function value.
- Otherwise, if the depth of recursion is less than some limit:
   - divide the pane into four quarters and invoke TransformPane
recursively on each quarter pane, gathering the returned wcs values
together and returning them to the caller.
- Otherwise:
   - Use the full mapping to transform the complete pane.




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

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> 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
> >>> http://mail.scipy.org/mailman/listinfo/astropy
> >
> >
> _______________________________________________
> AstroPy mailing list
> 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/4b772cb3/attachment.html>