Permalink
Browse files

Merge pull request #1421 from registerrier/map_mask

Add map region mask
  • Loading branch information...
registerrier committed Jun 3, 2018
2 parents d544802 + 078d533 commit 8ab962fa8b3e72e38b92318ec75d8ffaf4c28f5b
Showing with 72 additions and 1 deletion.
  1. +1 −0 gammapy/maps/base.py
  2. +3 −1 gammapy/maps/tests/test_geom.py
  3. +16 −0 gammapy/maps/tests/test_wcsnd.py
  4. +28 −0 gammapy/maps/wcs.py
  5. +24 −0 gammapy/maps/wcsnd.py
View
@@ -740,3 +740,4 @@ def set_by_idx(self, idx, vals):
Values vector. Pixels at `idx` will be set to these values.
"""
pass
@@ -4,7 +4,8 @@
from collections import OrderedDict
import numpy as np
from numpy.testing import assert_allclose
from astropy.coordinates import SkyCoord
from astropy.coordinates import SkyCoord, Angle
from regions import CircleSkyRegion
from ..geom import MapAxis, MapCoord
pytest.importorskip('scipy')
@@ -219,3 +220,4 @@ def test_mapcoords_to_coordsys():
'icrs').ra.deg, skycoord_gal.icrs.ra.deg)
assert_allclose(coords.skycoord.transform_to(
'icrs').dec.deg, skycoord_gal.icrs.dec.deg)
@@ -5,6 +5,7 @@
from numpy.testing import assert_allclose, assert_equal
from astropy.io import fits
from astropy.coordinates import SkyCoord
import astropy.units as u
from ...utils.testing import requires_dependency
from ..utils import fill_poisson
from ..geom import MapAxis, MapCoord, coordsys_to_frame
@@ -378,3 +379,18 @@ def test_coadd_unit():
assert_allclose(m1.data, 1.0001)
def test_make_region_mask():
from regions import CircleSkyRegion
geom = WcsGeom.create(npix=(3,3), binsz=2,
proj='CAR', coordsys='GAL')
m = WcsNDMap(geom)
region = CircleSkyRegion(SkyCoord(0, 0, unit='deg', frame='galactic'), 1.0*u.deg)
maskmap = m.make_region_mask(region)
assert maskmap.data.dtype == bool
assert np.sum(maskmap.data) == 1
maskmap = m.make_region_mask(region, inside=False)
assert np.sum(maskmap.data) == 8
View
@@ -8,6 +8,7 @@
from astropy.coordinates import SkyCoord
from astropy.coordinates.angle_utilities import angular_separation
import astropy.units as u
from regions import SkyRegion
from ..image.utils import make_header
from ..utils.scripts import make_path
from ..utils.wcs import get_resampled_wcs
@@ -634,6 +635,33 @@ def solid_angle(self):
return u.Quantity(dx * dy, 'sr')
def get_region_mask_array(self, region):
"""Return mask of pixels inside region in the the form of boolean array
Parameters
----------
region : `~regions.PixelRegion` or `~regions.SkyRegion` object
A region on the sky could be defined in pixel or sky coordinates.
Return
------
mask_array : `~numpy.ndarray` of booleans
the array of booleans
"""
from regions import PixCoord
# TODO : if Pixel Compound regions are taken into account, rather convert to PixelRegion
if isinstance(region, SkyRegion):
region = region.to_pixel(self.wcs)
if not self.is_regular:
raise NotImplementedError("Region mask is not working yet with multiresolution maps")
idx = self.get_idx()
pixcoord = PixCoord(idx[0],idx[1])
return region.contains(pixcoord)
def create_wcs(skydir, coordsys='CEL', projection='AIT',
cdelt=1.0, crpix=1., axes=None):
View
@@ -508,3 +508,27 @@ def plot(self, ax=None, idx=None, **kwargs):
im = ax.imshow(data, **kwargs)
ax.coords.grid(color='w', linestyle=':', linewidth=0.5)
return fig, ax, im
def make_region_mask(self, region, inside=True):
"""Create a mask of a given region
TODO: implement list of region for each axis
Parameters
----------
region : `~regions.PixelRegion` or `~regions.SkyRegion` object
A region on the sky could be defined in pixel or sky coordinates.
inside : bool
Output map is True inside the input region if inside is set to True and False outside and conversely.
Return
------
mask_map : `~gammapy.maps.WcsNDMap`
the mask map
"""
mask = self.geom.get_region_mask_array(region)
if inside is False:
np.logical_not(mask,out=mask)
# TODO : update meta table to include something about the region used for mask creation?
return WcsNDMap(geom=self.geom, data=mask, meta=self.meta)

0 comments on commit 8ab962f

Please sign in to comment.