<div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">On 17 October 2013 12:28, Maik Riechert <span dir="ltr"><<a href="mailto:maik.riechert@arcor.de" target="_blank">maik.riechert@arcor.de</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div bgcolor="#FFFFFF" text="#000000">
<br>
<div>On 17/10/13 08:54, David Berry wrote:<br>
</div>
<blockquote type="cite">
<div dir="ltr">AST attempts to speed up this sort of thing by
using a recursive scheme as follows:
<div class="gmail_extra">
<div class="gmail_quote">...
<div>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. </div>
<div><br>
</div>
<div>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.</div>
<div><br>
</div>
<div>See for instance <a href="http://starlink.jach.hawaii.edu/devdocs/sun211.htx/node430.html" target="_blank">http://starlink.jach.hawaii.edu/devdocs/sun211.htx/node430.html</a></div>
</div>
</div>
</div>
</blockquote>
<br>
I tried to use trangrid (with starlink-pyast) but I must be doing
something wrong:<br>
<br>
width = hdulist[0].header['NAXIS1']<br>
height = hdulist[0].header['NAXIS2']<br>
(wcs,_) = starlink.Atl.readfitswcs(hdulist[0])<br>
</div></blockquote><div><br></div><div style>Hi Maik,</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div bgcolor="#FFFFFF" text="#000000">mapping = wcs.getmapping( starlink.Ast.BASE, starlink.Ast.CURRENT )<br>
</div></blockquote><div><br></div><div style>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". </div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div bgcolor="#FFFFFF" text="#000000">
print mapping.trangrid([0,0],[width-1,height-1])<br>
<br>
[[-0.69877269 -0.69914315 -0.69951373 ..., -1.22724324 -1.22756251<br>
-1.22788175]<br>
[ 0.83534859 0.83540118 0.83545372 ..., 0.5841799 0.58413968<br>
0.58409942]]<br>
<br>
I am expecting coordinates ranging from 320..285 for RA and 48..33
for Dec.<br>
<br></div></blockquote><div><br></div><div style>Two issues:</div><div style><br></div><div style>1) AST uses radians to represent angles, not degrees. The above values are consistent with the ranges you quote taking this into account. </div>
<div style><br></div><div style>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).</div>
<div style><br></div><div style>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].</div>
<div style><br></div><div style>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).</div>
<div style><br></div><div style><br></div><div style>David <br></div><div style><br></div><div style><br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div bgcolor="#FFFFFF" text="#000000">
For comparison, my astropy code looks like that:<br>
<br>
wcs = WCS(hdulist[0].header)<br>
x, y = np.meshgrid(np.arange(width), np.arange(height)) <br>
ra, dec = wcs.wcs_pix2world(x, y, 0)<br>
<br>
Cheers<br>
Maik<br>
<br>
<blockquote type="cite">
<div dir="ltr">
<div class="gmail_extra">
<div class="gmail_quote">
<div>
David</div>
<div><br>
</div>
<div> </div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">Cheers,<br>
Tom<br>
<br>
><br>
> Maik<br>
><br>
><br>
> On 16/10/13 15:15, Thomas Robitaille wrote:<br>
>><br>
>> Hi Maik,<br>
>><br>
>> Someone asked a very similar question just
yesterday on GitHub:<br>
>><br>
>> <a href="https://github.com/astropy/astropy/issues/1587" target="_blank">https://github.com/astropy/astropy/issues/1587</a><br>
>><br>
>> In short, the solution in your case is:<br>
>><br>
>> NAXIS1 = 4256<br>
>> NAXIS2 = 2832<br>
>> x = np.arange(NAXIS1)<br>
>> y = np.arange(NAXIS2)<br>
>> X, Y = np.meshgrid(x, y)<br>
>> ra, dec = wcs.wcs_pix2world(X, Y, 0)<br>
>><br>
>> Any luck?<br>
>><br>
>> Tom<br>
>><br>
>><br>
>><br>
>><br>
>> On 16 October 2013 14:48, Maik Riechert <<a href="mailto:maik.riechert@arcor.de" target="_blank">maik.riechert@arcor.de</a>>
wrote:<br>
>>><br>
>>> Hi,<br>
>>><br>
>>> I'm quite new to python and astropy and got
stuck at the following which<br>
>>> probably is extremely easy to do.<br>
>>><br>
>>> I created a WCS object out of a FITS header
and wanted to get the RA/Dec<br>
>>> coordinates of each image pixel. I tried:<br>
>>><br>
>>> w.wcs_pix2world(np.ndindex(4256, 2832), 0)<br>
>>><br>
>>> where 4256x2832 are the image dimensions.
This returns 'ValueError:<br>
>>> object of too small depth for desired array'.
I was trying to avoid<br>
>>> allocating an array with all pixel
coordinates and thought an iterator<br>
>>> would work too.<br>
>>><br>
>>> What would you recommend for my case
(translating every pixel)?<br>
>>><br>
>>> Cheers and thanks,<br>
>>> Maik<br>
>>><br>
>>><br>
>>>
_______________________________________________<br>
>>> AstroPy mailing list<br>
>>> <a href="mailto:AstroPy@scipy.org" target="_blank">AstroPy@scipy.org</a><br>
>>> <a href="http://mail.scipy.org/mailman/listinfo/astropy" target="_blank">http://mail.scipy.org/mailman/listinfo/astropy</a><br>
><br>
><br>
_______________________________________________<br>
AstroPy mailing list<br>
<a href="mailto:AstroPy@scipy.org" target="_blank">AstroPy@scipy.org</a><br>
<a href="http://mail.scipy.org/mailman/listinfo/astropy" target="_blank">http://mail.scipy.org/mailman/listinfo/astropy</a><br>
</blockquote>
</div>
<br>
</div>
</div>
</blockquote>
<br>
</div>
</blockquote></div><br></div></div>