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