Add WCS keyword to provide WCS directly.

pyregion/core.py

@@ -45,14 +45,17 @@ def check_imagecoord(self):
         else:
             return True
 
-    def as_imagecoord(self, header):
+    def as_imagecoord(self, header, wcs=None):
         """New shape list in image coordinates.
 
         Parameters
         ----------
         header : `~astropy.io.fits.Header`
             FITS header
-
+
+        wcs : `~astropy.wcs.WCS`, optional
+            Full image WCS.  If not specified, will be parsed from `header`.
+
         Returns
         -------
         shape_list : `ShapeList`
@@ -66,7 +69,7 @@ def as_imagecoord(self, header):
             comment_list = cycle([None])
 
         r = RegionParser.sky_to_image(zip(self, comment_list),
-                                      header)
+                                      header, wcs=wcs)
         shape_list, comment_list = zip(*list(r))
         return ShapeList(shape_list, comment_list=comment_list)
 
@@ -88,7 +91,7 @@ def get_mpl_patches_texts(self, properties_func=None,
 
         return patches, txts
 
-    def get_filter(self, header=None, origin=1):
+    def get_filter(self, header=None, origin=1, wcs=None):
         """Get filter.
         Often, the regions files implicitly assume the lower-left
         corner of the image as a coordinate (1,1). However, the python
@@ -103,7 +106,10 @@ def get_filter(self, header=None, origin=1):
             FITS header
         origin : {0, 1}
             Pixel coordinate origin
-
+        wcs : `~astropy.wcs.WCS`, optional
+            Full image WCS.  If not specified, will be parsed from `header if
+            necessary.
+
         Returns
         -------
         filter : TODO
@@ -117,13 +123,13 @@ def get_filter(self, header=None, origin=1):
                 raise RuntimeError("the region has non-image coordinate. header is required.")
             reg_in_imagecoord = self
         else:
-            reg_in_imagecoord = self.as_imagecoord(header)
+            reg_in_imagecoord = self.as_imagecoord(header, wcs=wcs)
 
         region_filter = as_region_filter(reg_in_imagecoord, origin=origin)
 
         return region_filter
 
-    def get_mask(self, hdu=None, header=None, shape=None):
+    def get_mask(self, hdu=None, header=None, shape=None, wcs=None):
         """Create a 2-d mask.
 
         Parameters
@@ -134,6 +140,9 @@ def get_mask(self, hdu=None, header=None, shape=None):
             FITS header
         shape : tuple
             Image shape
+        wcs : `~astropy.wcs.WCS`, optional
+            Full image WCS.  If not specified, will be parsed from `header` if 
+            necessary.
 
         Returns
         -------
@@ -152,7 +161,7 @@ def get_mask(self, hdu=None, header=None, shape=None):
         if hdu and shape is None:
             shape = hdu.data.shape
 
-        region_filter = self.get_filter(header=header)
+        region_filter = self.get_filter(header=header, wcs=wcs)
         mask = region_filter.mask(shape)
 
         return mask

pyregion/ds9_region_parser.py

@@ -159,7 +159,7 @@ def convert_attr(self, l):
                 yield l1, c1
 
     @staticmethod
-    def sky_to_image(shape_list, header):
+    def sky_to_image(shape_list, header, wcs=None):
         """Converts a `ShapeList` into shapes with coordinates in image coordinates
 
         Parameters
@@ -168,7 +168,10 @@ def sky_to_image(shape_list, header):
             The ShapeList to convert
         header : `~astropy.io.fits.Header`
             Specifies what WCS transformations to use.
-
+        wcs : `~astropy.wcs.WCS`, optional
+            Full image WCS.  If not specified, will be parsed from `header` if 
+            necessary.
+
         Yields
         -------
         shape, comment : Shape, str
@@ -184,7 +187,7 @@ def sky_to_image(shape_list, header):
             if isinstance(shape, Shape) and \
                     (shape.coord_format not in image_like_coordformats):
 
-                new_coords = convert_to_imagecoord(shape, header)
+                new_coords = convert_to_imagecoord(shape, header, wcs=wcs)
 
                 l1n = copy.copy(shape)
 

pyregion/tests/test_get_mask.py

@@ -1,7 +1,9 @@
 import os
 import numpy as np
 from os.path import join
-from astropy.io.fits import Header
+from astropy.io.fits import Header, ImageHDU
+import astropy.wcs
+
 from .. import open as pyregion_open
 
 rootdir = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'data')
@@ -21,3 +23,19 @@ def test_region():
     assert isinstance(mask, np.ndarray) and mask.shape == (100, 100)
 
     # TODO: assert the content of the mask, too
+
+def test_region_wcs():
+    ref_name = "test01_fk5.reg"
+
+    header = demo_header()
+
+    ref_region = pyregion_open(join(rootdir, ref_name)).as_imagecoord(header)
+
+    wcs = astropy.wcs.WCS(header)
+
+    shape = (header['NAXIS2'], header['NAXIS1'])
+
+    mask_from_header = ref_region.get_mask(header=header, shape=shape)
+    mask_from_wcs = ref_region.get_mask(shape=shape, wcs=wcs)
+
+    assert (mask_from_header.sum() == mask_from_wcs.sum()) and (mask_from_header.shape == shape)

pyregion/tests/test_region.py

@@ -3,6 +3,8 @@
 import numpy as np
 from os.path import join
 from astropy.io.fits import Header
+import astropy.wcs
+
 from .. import open as pyregion_open
 from numpy.testing import assert_allclose
 
@@ -55,6 +57,50 @@ def test_region(ref_name, reg_name, header_name):
 
         assert ref_reg.exclude == reg.exclude
 
+@pytest.mark.parametrize(("ref_name", "reg_name", "header_name"), [
+    ("test01_img.reg", "test01_fk5_sexagecimal.reg", "sample_fits01.header"),
+    ("test01_img.reg", "test01_gal.reg", "sample_fits01.header"),
+    ("test01_img.reg", "test01_ds9_physical.reg", "sample_fits01.header"),
+    ("test01_img.reg", "test01_fk5_degree.reg", "sample_fits01.header"),
+    ("test01_img.reg", "test01_mixed.reg", "sample_fits01.header"),
+    ("test01_img.reg", "test01_ciao.reg", "sample_fits01.header"),
+    ("test01_img.reg", "test01_ciao_physical.reg", "sample_fits01.header"),
+    ("test01_img.reg", "test01_fk5.reg", "sample_fits01.header"),
+    ("test01_img.reg", "test01_fk4.reg", "sample_fits01.header"),
+    ("test01_img.reg", "test01_icrs.reg", "sample_fits01.header"),
+    ("test02_1_img.reg", "test02_1_fk5.reg", "sample_fits02.header"),
+    ("test_annuli.reg", "test_annuli_wcs.reg", "sample_fits01.header"),
+    ("test03_img.reg", "test03_fk5.reg", "sample_fits03.header"),
+    ("test03_img.reg", "test03_icrs.reg", "sample_fits03.header"),
+    ("test03_img.reg", "test03_ciao_physical.reg", "sample_fits03.header"),
+    ("test03_img.reg", "test03_gal.reg", "sample_fits03.header"),
+])
+def test_region_wcs(ref_name, reg_name, header_name):
+    header = Header.fromtextfile(join(rootdir, header_name))
+    wcs = astropy.wcs.WCS(header=header)
+
+    # Header-only
+    ref_region = pyregion_open(join(rootdir, ref_name)).as_imagecoord(header)
+
+    # With WCS
+    r = pyregion_open(join(rootdir, reg_name)).as_imagecoord(header, wcs=wcs)
+
+    assert len(r) == len(ref_region)
+
+    for ref_reg, reg in zip(ref_region, r):
+        if reg.name == "rotbox":
+            reg.name = "box"
+
+        assert ref_reg.name == reg.name
+
+        # Normalize everything like angles
+        ref_list = np.asarray(ref_reg.coord_list)
+        reg_list = np.asarray(reg.coord_list)
+        assert_allclose((ref_list + 180) % 360 - 180,
+                        (reg_list + 180) % 360 - 180,
+                        atol=0.03)
+
+        assert ref_reg.exclude == reg.exclude
 
 @pytest.mark.parametrize("reg_name", [
     "test_annuli_ciao.reg",  # subset of test03_img.reg

pyregion/wcs_converter.py

@@ -55,7 +55,7 @@ def _generate_arg_types(coordlist_length, shape_name):
     return arg_types
 
 
-def convert_to_imagecoord(shape, header):
+def convert_to_imagecoord(shape, header, wcs=None):
     """Convert the coordlist of `shape` to image coordinates
 
     Parameters
@@ -66,6 +66,12 @@ def convert_to_imagecoord(shape, header):
     header : `~astropy.io.fits.Header`
         Specifies what WCS transformations to use.
 
+    wcs : `~astropy.wcs.WCS`, optional
+        Full image WCS.  If not specified, will be parsed from `header` if 
+        necessary, e.g., with
+
+           >>> new_wcs = astropy.wcs.WCS(header)
+
     Returns
     -------
     new_coordlist : list
@@ -78,7 +84,11 @@ def convert_to_imagecoord(shape, header):
     is_even_distance = True
     coord_list_iter = iter(zip(shape.coord_list, arg_types))
 
-    new_wcs = WCS(header)
+    if wcs is None:
+        new_wcs = WCS(header)
+    else:
+        new_wcs = wcs
+
     pixel_scales = proj_plane_pixel_scales(new_wcs)
 
     for coordinate, coordinate_type in coord_list_iter: