[AstroPy] Running wcs_pix2world for all image pixels

Thu Oct 17 08:41:19 EDT 2013

On 17 October 2013 12:28, Maik Riechert <maik.riechert at arcor.de> wrote:

>
> 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])
>

Hi Maik,

> mapping = wcs.getmapping( starlink.Ast.BASE, starlink.Ast.CURRENT )
>

You don't really need this invocation of getmapping. The "wcs" object is a
FrameSet, which is a subclass of Mapping. When used as a Mapping, a
FrameSet object behaves like the Mapping from its base Frame to its current
frame, so you can apply the trangrid method directly to "wcs" instead of
"mapping".

> 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.
>
>
Two issues:

1) AST uses radians to represent angles, not degrees. The above values are
consistent with the ranges you quote taking this into account.

2) The negative RA values need 2*PI added to them. Longitude axes such as
RA have the inconvenient characteristic of having a wrap-around
discontinuity somewhere. When reporting values to human beings is is usual
to put this discontinuity at RA=0. But this is not always what you want
when dealing with RA values in code. For instance, if your image contains
zero longitude, measuring RA intervals on your image is awkward. In this
case it would be much more convenient to put the discontinuity at RA=180,
so that RA=-10 is used in place of RA=350 (for instance).

In view of such uncertainty about the best place to put the discontinuity,
AST takes a neutral view and makes no assumption or guarantee about where
the discontinuity will be. Instead, it provides the astNorm method (a
method of the Frame class)  which "normalises" the axis values representing
a position within the Frame, to the form usually expected by human readers.
So for instance the SkyFrame class provides an implementation of astNorm
that puts longitude into the range [0,2*PI] and latitude into the range
[-PI/2,+PI/2].

AST's expectation is that the bulk of calculations performed on
(longitude,latitude) pairs will be unaffected by the position of the
discontinuity, and so it saves time by not performing normalisation
automatically. In the few cases where it really does matter that the
discontinuity is at a specific place,  the astNorm method can be used (or
you can do the normalisation using your own code if you know a priori that
the axis values represents a position on the sky).

David

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>
>> 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/df4ac3fa/attachment.html>