Spatial filter seems to match positions outside region

[grahambell]

The spatial filter seems to match positions which are not part of the specified region. In the following (very slow) test script, I test points over the whole sky using a 1 degree grid. The output has, as well as the two expected circles, lines at RA = +/- 90 degrees and a few other scattered points. This was with Astropy 1.0.6 and the current master version of pyregion from this repository.

from __future__ import print_function                                           
from astropy.units import degree                                                
from astropy.wcs import WCS                                                     
import matplotlib.pyplot as plt                                                 
import pyregion                                                                 

region = pyregion.parse('''                                                     
    fk5                                                                         
    circle(0:00:00.000,00:00:00.00,36000")                                      
    circle(3:00:00.000,75:00:00.00,36000")                                      
''')                                                                            

ras = []                                                                        
decs = []                                                                       

for ra in range(-179, 180):                                                     
    for dec in range (-89, 90):                                                 
        wcs = WCS(naxis=2)                                                      
        wcs.wcs.radesys = 'ICRS'                                                
        wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN']                                
        wcs.wcs.cunit = [degree, degree]                                        
        wcs.wcs.crpix = [1, 1]                                                  
        wcs.wcs.cdelt = [-0.0003, 0.0003]                                       
        wcs.wcs.crval = [ra, dec]                                               

        if region.get_filter(header=wcs.to_header()).inside1(1, 1):             
            #print(ra, dec)                                                     
            ras.append(ra)                                                      
            decs.append(dec)                                                    

plt.scatter(ras, decs)                                                          
plt.xlim(-180, 180)                                                             
plt.ylim(-90, 90)                                                               
plt.show() 

[cdeil]

@grahambell - Thanks for the bug report. I see the same output with Python 3.5, Astropy 1.2, Numpy 1.11.2, Cython 0.24.1.

@grahambell - By any chance, do you already have a test case for this that can be run quickly? E.g. same example, but just with one or a few RA / DEC values where the result is incorrect?
If I don't hear back I'll start investigating this after lunch.

For now I've put a v1.3 milestone on this.
@leejjoon - If you consider this a blocker for the v1.2 release this week, please change the milestone to v1.2.

[cdeil]

I've started debugging this ... here's some notes.

Script to print the RA / DEC values where the pyregion result is incorrect:
https://gist.github.com/cdeil/ee69188cb8df95a4b8e2fca22dceea5c

First incorrect result is for:

RA, DEC = -141 15

Script to show what's going on for that case:
https://gist.github.com/cdeil/91ef8dedaa9e9fafc4e460ba1355c6af

Output:
https://gist.github.com/cdeil/5236392a1375006e98af55a2a198211f

Looks like the issue is nan values in the first circle in the filter = region.get_filter(header=header):

Or[Circle(nan, nan, nan), Circle(3772783.136359, 139402848.023292, 147909026.067737)]

However, this seems to work: region.as_imagecoord(header=header)

[Shape : circle ( HMS(0:00:00.000),DMS(00:00:00.00),Ang(36000") ), Shape : circle ( HMS(3:00:00.000),DMS(75:00:00.00),Ang(36000") )]

I'll continue looking at this after lunch.

@leejjoon or anyone -- any idea what the issue or fix is?

[cdeil]

Changing milestone for this issue to 1.2.

[cdeil]

I've narrowed this down to a bug (I think) in RegionParser.sky_to_image.

Here's a small self-contained test case:

from astropy.wcs import WCS
wcs = WCS(naxis=2)
wcs.wcs.radesys = 'ICRS'
wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN']
wcs.wcs.cunit = ['deg', 'deg']
wcs.wcs.crpix = [1, 1]
wcs.wcs.cdelt = [-0.0003, 0.0003]
wcs.wcs.crval = [-141, 15]
header = wcs.to_header()


import pyregion
shape_list = pyregion.parse('fk5\ncircle(0:00:00.000,00:00:00.00,36000")')
# import IPython; IPython.embed()
shape_list2 = shape_list.as_imagecoord(header=header)

print(repr(header))
print('shape_list: ', shape_list)
print('shape_list2:', shape_list2)

Output:

WCSAXES =                    2 / Number of coordinate axes                      
CRPIX1  =                  1.0 / Pixel coordinate of reference point            
CRPIX2  =                  1.0 / Pixel coordinate of reference point            
CDELT1  =              -0.0003 / [deg] Coordinate increment at reference point  
CDELT2  =               0.0003 / [deg] Coordinate increment at reference point  
CUNIT1  = 'deg'                / Units of coordinate increment and value        
CUNIT2  = 'deg'                / Units of coordinate increment and value        
CTYPE1  = 'RA---TAN'           / Right ascension, gnomonic projection           
CTYPE2  = 'DEC--TAN'           / Declination, gnomonic projection               
CRVAL1  =               -141.0 / [deg] Coordinate value at reference point      
CRVAL2  =                 15.0 / [deg] Coordinate value at reference point      
LONPOLE =                180.0 / [deg] Native longitude of celestial pole       
LATPOLE =                 15.0 / [deg] Native latitude of celestial pole        
RADESYS = 'ICRS'               / Equatorial coordinate system                   
shape_list:  [Shape : circle ( HMS(0:00:00.000),DMS(00:00:00.00),Ang(36000") )]
shape_list2: [Shape : circle ( HMS(0:00:00.000),DMS(00:00:00.00),Ang(36000") )]

So for this header and shape_list, ShapeList.as_imagecoord doesn't convert to image coordinates or raise an error, but just return the shape in sky coordinates again.

I've looked at bit at the RegionParser.sky_to_image implementation. It executes the first branch of the if statement here. But I don't know what's going on ... that's a very complicated 100 lines long method that's hard to test because it does so much in one function.

I'll set the milestone back to 1.3 ... probably that bug has been there for years and IMO this shouldn't hold up the 1.2 release, which is otherwise ready to go.

@leejjoon or anyone else: do you have time to look into this issue?

[cdeil]

Here's the full list of pixels (RA, DEC) in the first two columns where pyregion currently gives incorrect results for the original test script (takes ~ 20 minutes to run, so I'm pasting it here for reference):
https://gist.github.com/cdeil/b38e3f571ebfe738582634d8243338f5

[leejjoon]

I try to track this problem, and the script below reproduces the source of the problem.

from astropy.wcs import WCS
from astropy.coordinates import SkyCoord

wcs = WCS(naxis=2)
wcs.wcs.radesys = 'ICRS'
wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN']
wcs.wcs.cunit = ['deg', 'deg']
wcs.wcs.crpix = [1, 1]
wcs.wcs.cdelt = [-0.0003, 0.0003]
wcs.wcs.crval = [-141, 15]
header = wcs.to_header()
new_wcs = WCS(header)

c1, c2 = 0, 0
coord_format = "fk5"
old_coordinate = SkyCoord(c1, c2,
                          frame=coord_format, unit='degree',
                          obstime='J2000')

print(old_coordinate.to_pixel(new_wcs, origin=1))

The resulting coordinates are NaNs.

As far as I can see, this is because the coordinate given ([0, 0] in this example) is either not defined in the given projection or it overflows.

I note that region filters that pyregion generates are based on their image coordinate. And things like radius is converted to a radius in the image plane simply using the plate scale. So, you should not use pyregion's filter if the sky curvature is not negligible.

For pyregion, I guess the question is what should be its behavior when some of the input points are not defined in the image coordinate. My inclination is that we just raise an exception indicating that the filter is ill-defined. Thought?

[cdeil]

Maybe @nden or @astrofrog or someone that knows astropy.wcs well can comment:
In the example given by @leejjoon in the comment above this one, is the nan output expected / correct?

[cdeil]

Is this a blocker for a pyregion 2.0 release or can it be moved to 2.1?

@astrofrog or @keflavich - Could you please have a look at the example by @leejjoon above and comment whether this is expected behaviour for this WCS or a bug?

For pyregion, I guess the question is what should be its behavior when some of the input points are not defined in the image coordinate. My inclination is that we just raise an exception indicating that the filter is ill-defined. Thought?

I think for operations that act on arrays the common idiom is to return NaN for items where no result is available. That's e.g. what you get from WCS for pixels where the transformation isn't defined. So I'm -1 on raising an exception, because that means users don't get any output. Numpy emits warnings:

>>> import numpy as np
>>> np.sqrt([3, -2])
RuntimeWarning: invalid value encountered in sqrt
array([ 1.73205081,         nan])

which in my experience are more often annoying and have to be suppressed instead of helpful, but emitting warnings is certainly also a valid way to deal with this.

[keflavich]

If the projection diverges, e.g., if you try to use a tangent projection at one of the poles, I think it would be nice to raise a helpful exception. This is a case where an answer is theoretically possible, but it is definitely wrong, whereas the case @cdeil highlighted is where there is no theoretical answer; I'm happy with a NaN return in that situation.

However, looking at the specific case @leejjoon posted and the original bug report, I don't know what's going on, since the requested coordinates that evaluate to NaN are generally nowhere near a pole. Or, does the TAN projection use the local coordinate as the origin? That would make sense... I'd have to do some more reading and/or testing to understand really what's going on.

[leejjoon]

I think for operations that act on arrays the common idiom is to return NaN for items where no result is available. That's e.g. what you get from WCS for pixels where the transformation isn't defined. So I'm -1 on raising an exception, because that means users don't get any output. Numpy emits warnings:

Then, how about changing the behavior of the filter so that it raises/warns if the input regions has NaN.

[leejjoon]

Or, does the TAN projection use the local coordinate as the origin?

Yes.

[leejjoon]

Sorry, closed accidentally. Reopening it.

[cdeil]

I also think that returning NaN is fine / preferable to warnings. As far as I know, this is consistent with the behaviour of astropy.wcs which simply returns NaN for pix_to_world where the pixel isn't on the sphere (such as e.g. at the edges for all-sky images in AIT projection).

As far as I can tell though, the spatial filter results for the TAN examples given above indicate a bug, i.e. simply aren't correct for unknown reasons (see image posted above). Is that true?

Is there anyone that has the time and expertise to have a look and resolve this issue?
(at the moment this is under the v2.0 milestone, but maybe we should just do the release now and move this to the next milestone?)

[cdeil]

The last stable pyregion 1.2 release is broken, so I'll go ahead and make a new release v2.0 now (see #94 )
Moving this issue to the v2.1 milestone now.