[AstroPy] Bug in wcs_world2pix?

Michael Droettboom mdroe at stsci.edu
Sat Sep 20 17:58:42 EDT 2014


With the help of Mark Calabretta convincing me those input RA and Dec 
should be rejected by bounds checking and marked as invalid, I 
discovered a long standing and strangely-never-noticed bug in astropy.  
It was checking for the wrong status code that indicates when one of the 
input world coordinates is invalid -- the status codes are different for 
wcsp2s and wcss2p, and I think I just didn't notice that.

A fix is available in the linked PR below, which should make it into the 
next bugfix release of astropy.  With this fix, your example will now 
return the pixel coordinates as "nan", which is how astropy.wcs 
indicates invalid values.

https://github.com/astropy/astropy/pull/2965

Mike

On 09/20/2014 07:23 AM, Maik Riechert wrote:
> I still think it's a bug, but anyway, I worked around it be converting
> to RA,Dec again. This is way easier than somehow calculating the
> distance for arbitrary WCS.
>
> The code is here for interested people:
> https://gist.github.com/neothemachine/2fd64d7a552ba7743ce7
>
> It's not directly usable (e.g. misses constellation data) but the code
> should give useful insights.
>
> An example image:
> http://oi62.tinypic.com/sncjkp.jpg
>
> (Note that I used padding of 10px around stars and clipped the drawing
> below the horizon.)
>
> Cheers and thanks!
> Maik
>
> Am 19.09.2014 17:22, schrieb Eduardo Gonzalez-Solares:
>> As others have said you cannot represent the whole sky using only one TAN projection. More than a bug it is a "feature". What I would do is check first that your constellation is at a roughly reasonable distance (depends on the size of your image) of the centre of your image. That would probably make your routine quicker too.
>>
>> Cheers,
>>
>> Eduardo
>>
>> On 19 Sep 2014, at 15:57, Maik Riechert <maik.riechert at arcor.de> wrote:
>>
>>> Just to add some context of what I'm doing. I got the star constellation
>>> data from Marshall Perrin the other day on this list. And now I'm trying
>>> to draw these constellations onto an image using WCS data. That means,
>>> for each constellation, I convert the RA,DEC points to pixel coordinates
>>> using wcs_world2pix, then I check whether all points are outside the
>>> image bounds (in that case I skip drawing it), if they are partly or all
>>> inside, then I draw the constellation.
>>>
>>> So yes, I think the best way would be to return NaN in these cases.
>>> Otherwise I wouldn't know how to easily filter out those wrong points.
>>>
>>> Cheers
>>> Maik
>>>
>>>
>>> On 19/09/14 15:14, Thomas Robitaille wrote:
>>>> Ah of course, this is indeed what is going on here! I agree that it's
>>>> worth thinking about whether cases like this should return an error.
>>>> In the past (in WCSAxes) I've simply checked if the coordinate
>>>> round-trips in order to determine if it's in the observational plane,
>>>> but it may be more sensible to simply return NaN.
>>>>
>>>> Cheers,
>>>> Tom
>>>>
>>>> On 19 September 2014 14:14, David Berry <d.berry at jach.hawaii.edu> wrote:
>>>>> Am I missing something? This is a TAN projection centred on CRVAL1,
>>>>> CRVAL2#. Since it is a TAN projection, the projection plane only
>>>>> includes half the sky - i.e. the hemisphere centred on (16,23). But
>>>>> the test point (211,-26) is not contained within this hemisphere and
>>>>> so is not included in the TAN projection plane. In other words, no
>>>>> pixel has the coordinates RA!1 Dec=6, and so the results of  using
>>>>> wcs.wcs_world2pix should be NaN, or "invalid", or something.
>>>>>
>>>>> Unless I'm missing something...
>>>>>
>>>>> David
>>>>>
>>>>>
>>>>> On 19 September 2014 12:58, Bob Garwood <bgarwood at nrao.edu> wrote:
>>>>>> I suspect what's going on here has to do with the specific CD matrix, which
>>>>>> includes rotation. I think because of that rotation the world coordinate
>>>>>> actually maps back in to the grid at the antipode. I haven't checked this by
>>>>>> hand to see if that's the case.
>>>>>>
>>>>>> Bob
>>>>>>
>>>>>> On September 19, 2014 7:47:30 AM EDT, Howard Bushouse <bushouse at stsci.edu>
>>>>>> wrote:
>>>>>>> Huh. But that’s not how astronomers measure RA and Dec, even when
>>>>>>> translated to degrees. If I’ve got an object with an RA of 17:45:00.0
>>>>>>> (hh:mm:ss.w) the RA in degrees is 266.25. We *never* use a negative RA.
>>>>>>> Instead of a range of -180 to 180, we use 0 to 360. Dec is the only one that
>>>>>>> ever goes negative (-90 to +90, as you indicated).
>>>>>>>
>>>>>>> -hb
>>>>>>>
>>>>>>> On Sep 18, 2014, at 8:26 PM, Michael Droettboom <mdroe at stsci.edu> wrote:
>>>>>>>
>>>>>>>>   The range of ra and dec is (-180, 180) and (-90, 90) respectively.
>>>>>>>>   (211, -26) is equivalent to (31, 26).
>>>>>>>>
>>>>>>>>   Mike
>>>>>>>>
>>>>>>>>   On 09/18/2014 05:04 PM, Maik Riechert wrote:
>>>>>>>>>   Hi,
>>>>>>>>>
>>>>>>>>>   I'm having some great trouble currently and am a bit confused to!
>>>>>>>>>   o:
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>   wcs =(header)
>>>>>>>>>   x,y =.wcs_world2pix(211, -26, 0)
>>>>>>>>>   print x,y # 836.316942718 658.26364248
>>>>>>>>>   ra,dec =.wcs_pix2world(x, y, 0)
>>>>>>>>>   print ra, dec # 31.0 26.0
>>>>>>>>>
>>>>>>>>>   It seems like world2pix gives me bogus results. I was expecting values
>>>>>>>>>   way outside my image area (image size is 4256x2832). Am I doing
>>>>>>>>>   something obvious wrong?
>>>>>>>>>
>>>>>>>>>   The WCS headers can be seen at http://pastebin.com/JrdiL949
>>>>>>>>>
>>>>>>>>>   Thanks for the help,
>>>>>>>>>   Maik
>>>>>>>>> ________________________________
>>>>>>>>>
>>>>>>>>>   AstroPy mailing list
>>>>>>>>>   AstroPy at scipy.org
>>>>>>>>>   http://mail.scipy.org/mailman/listinfo/astropy
>>>>>>>>
>>>>>>>>   --
>>>>>>>>   Michael Droettboom
>>>>>>>>   Science Software Branch
>>>>>>>>   Space Telescope Science Institute
>>>>>>>>
>>>>>>>>   http://www.droettboom.com
>>>>>>>>
>>>>>>>> ________________________________
>>>>>>>>
>>>>>>>>   AstroPy mailing list
>>>>>>>>   AstroPy at scipy.or!
>>>>>>>>   g
>>>>>>>>
>>>>>>>> http://mail.scipy.org/mailman/listinfo/astropy
>>>>>>> ________________________________
>>>>>>>
>>>>>>> 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
>>>>>>
>>>>> _______________________________________________
>>>>> 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
>>
> _______________________________________________
> AstroPy mailing list
> AstroPy at scipy.org
> http://mail.scipy.org/mailman/listinfo/astropy


-- 
Michael Droettboom
Science Software Branch
Space Telescope Science Institute

http://www.droettboom.com




More information about the AstroPy mailing list