[AstroPy] Running wcs_pix2world for all image pixels
Maik Riechert
maik.riechert at arcor.de
Thu Oct 17 11:28:45 EDT 2013
On 17/10/13 17:20, David Berry wrote:
>
>
>
> On 17 October 2013 15:28, Maik Riechert <maik.riechert at arcor.de
> <mailto:maik.riechert at arcor.de>> wrote:
>
>
> On 17/10/13 14:41, David Berry wrote:
>> On 17 October 2013 12:28, Maik Riechert <maik.riechert at arcor.de
>> <mailto:maik.riechert at arcor.de>> wrote:
>>
>>
>> I tried to use trangrid (with starlink-pyast) but I must be
>> doing something wrong:
>> ...
>>
>> [[-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).
>>
>
> Thanks David for this detailed explanation. I could confirm the
> results now too. When using for instance 0.00005 as tolerance it
> is around 10 times faster than the exact algorithm, very nice!
>
> One more question though. I translated three single pixel
> coordinates to RA, Dec and then again to pixel coordinates using
> the tran() method. I expected that both sets of pixel coordinates
> are nearly identical (which they are when using astropy), but this
> is what I got:
>
> Input pixel coords:
> [[ -196., 735., 4059.],
> [ -504., 459., 2328.]]
>
> Calculated pixel values:
> [[ -205.93229216, 735. , 4179.93253144],
> [ -499.36855108, 459. , 2246.86209704]]
>
> Note that (735,459) is the WCS center reference.
>
> Any idea?
>
>
>
> Hi Maik,
> The tolerance you specify when running TranGrid is a
> distance in the output coordinate system, which in your case I presume
> is (ra,dec). So a tolerance value of 0.00005 is taken as a value in
> radians, which is about 10 arc-seconds. So errors of up to this sort
> of size could be expected from the linear approximations being used.
> What is the pixel size in your case? You should find that using a
> smaller tolerance results in smaller round-trip errors. If this is not
> the case, it may be a bug, in which case I'd be interested to see the
> FITS header you are using so that I can track it down.
I didn't use tranGrid here, just tran! So no tolerance involved.
NAXIS1 = 1468 / image width
NAXIS2 = 916 / image height
CTYPE1 = 'RA---TAN-SIP' / TAN (gnomic) projection + SIP distortions
CTYPE2 = 'DEC--TAN-SIP' / TAN (gnomic) projection + SIP distortions
EQUINOX = 2000.0 / Equatorial coordinates definition (yr)
LONPOLE = 180.0 / no comment
LATPOLE = 0.0 / no comment
CRVAL1 = 302.88771538 / RA of reference point
CRVAL2 = 41.8941029213 / DEC of reference point
CRPIX1 = 735 / X reference pixel
CRPIX2 = 459 / Y reference pixel
CUNIT1 = 'deg ' / X pixel scale units
CUNIT2 = 'deg ' / Y pixel scale units
CD1_1 = -0.0157184611387 / Transformation matrix
CD1_2 = 0.00022032143241 / no comment
CD2_1 = -0.000284856284919 / no comment
CD2_2 = -0.0159943489114 / no comment
IMAGEW = 1468 / Image width, in pixels.
IMAGEH = 916 / Image height, in pixels.
A_ORDER = 2 / Polynomial order, axis 1
A_0_2 = -5.10167355829E-05 / no comment
A_1_1 = -5.54267287337E-05 / no comment
A_2_0 = 3.68095680847E-05 / no comment
B_ORDER = 2 / Polynomial order, axis 2
B_0_2 = 2.49853669312E-05 / no comment
B_1_1 = 6.07756179822E-05 / no comment
B_2_0 = -1.70350503136E-05 / no comment
AP_ORDER= 2 / Inv polynomial order, axis 1
AP_0_1 = -0.0014471978634 / no comment
AP_0_2 = 5.0935282476E-05 / no comment
AP_1_0 = 0.000641040653934 / no comment
AP_1_1 = 5.54267224252E-05 / no comment
AP_2_0 = -3.674588665E-05 / no comment
BP_ORDER= 2 / Inv polynomial order, axis 2
BP_0_1 = 0.0010860421788 / no comment
BP_0_2 = -2.49999638422E-05 / no comment
BP_1_0 = -0.000569392740582 / no comment
BP_1_1 = -6.07473585161E-05 / no comment
BP_2_0 = 1.7000827442E-05 / no comment
Maik
>
> David
>
>
> Maik
>
>
>>
>> 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 <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/64e950ed/attachment.html>
More information about the AstroPy
mailing list