diff --git a/pyregion/core.py b/pyregion/core.py index 5890018..11aed3f 100644 --- a/pyregion/core.py +++ b/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 diff --git a/pyregion/ds9_region_parser.py b/pyregion/ds9_region_parser.py index 1b1e82b..e3767a4 100644 --- a/pyregion/ds9_region_parser.py +++ b/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) diff --git a/pyregion/tests/test_get_mask.py b/pyregion/tests/test_get_mask.py index 95a6c7d..ace76c5 100644 --- a/pyregion/tests/test_get_mask.py +++ b/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) diff --git a/pyregion/tests/test_region.py b/pyregion/tests/test_region.py index 55c0789..a6f5a86 100644 --- a/pyregion/tests/test_region.py +++ b/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 diff --git a/pyregion/wcs_converter.py b/pyregion/wcs_converter.py index b9481b8..18b119f 100644 --- a/pyregion/wcs_converter.py +++ b/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: