test_issue_65.py

from __future__ import print_function
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
import matplotlib.pyplot as plt
import regions
import pyregion

region_string = '''
fk5
circle(0:00:00.000,00:00:00.00,36000")
circle(3:00:00.000,75:00:00.00,36000")
'''

region = pyregion.parse(region_string)

region2 = regions.ds9_string_to_objects(region_string)
region2 = region2[0] & region2[1]

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 = [u.degree, u.degree]
        wcs.wcs.crpix = [1, 1]
        wcs.wcs.cdelt = [-0.0003, 0.0003]
        wcs.wcs.crval = [ra, dec]

        inside_pyregion = region.get_filter(header=wcs.to_header()).inside1(1, 1)
        inside_regions = region2.contains(SkyCoord(ra, dec, frame='fk5', unit='deg'))

        if inside_pyregion != inside_regions:
            print(ra, dec, inside_pyregion, inside_regions)

        if inside_pyregion:
            ras.append(ra)
            decs.append(dec)

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