New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add map region mask #1421

Merged
merged 7 commits into from Jun 3, 2018
View
@@ -740,3 +740,4 @@ def set_by_idx(self, idx, vals):
Values vector. Pixels at `idx` will be set to these values.
"""
pass
View
@@ -867,6 +867,29 @@ def apply_mask(self, mask):
return self.__class__(data, self.coordsys,
match_by_name=self._match_by_name)
def apply_region_mask(self, region, inside=True):
"""Return a copy of this coordinate object after applying a region selection mask.
Warning: Selection is performed on spatial coordinates only, i.e. the same region is used
for all non-spatial axes.
Parameters
----------
region : `~region.SkyRegion`
Region defining the mask
inside : boolean
If True keep only coordinates inside, if False only coordinates outside
Returns
-------
coords : `~MapCoord`
A coordinates object.
"""
mask = region.contains(self.skycoord)
if inside is False:
np.logical_not(mask, out=mask)
return self.apply_mask(mask)
class MapGeomMeta(InheritDocstrings, abc.ABCMeta):
pass
View
@@ -10,6 +10,7 @@
from astropy.io import fits
from astropy.coordinates import SkyCoord
from astropy.units import Quantity
from regions import SkyRegion

This comment has been minimized.

@cdeil

cdeil May 29, 2018

Member

Remove unused import

@cdeil

cdeil May 29, 2018

Member

Remove unused import

This comment has been minimized.

@registerrier

registerrier Jun 1, 2018

Contributor

Done

@registerrier

registerrier Jun 1, 2018

Contributor

Done

from ..utils.scripts import make_path
from .wcs import WcsGeom
from .geom import MapGeom, MapCoord, pix_tuple_to_idx
View
@@ -4,6 +4,7 @@
import numpy as np
from astropy.io import fits
from astropy.units import Quantity
from regions import SkyRegion

This comment has been minimized.

@cdeil

cdeil May 29, 2018

Member

Remove unused import

@cdeil

cdeil May 29, 2018

Member

Remove unused import

This comment has been minimized.

@registerrier

registerrier Jun 1, 2018

Contributor

Done

@registerrier

registerrier Jun 1, 2018

Contributor

Done

from .utils import unpack_seq
from .geom import MapCoord, pix_tuple_to_idx, coord_to_idx
from .utils import interp_to_order
@@ -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,19 @@ 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)
def test_mapcoord_region_filter():
lon, lat = np.array([0.0, 10.0]), np.array([0.0, 10.0])
energy = np.array([100., 1000.])
coords = MapCoord.create(
dict(lon=lon, lat=lat, energy=energy), coordsys='GAL')
region = CircleSkyRegion(SkyCoord(0.,0.,frame='galactic',unit='deg'), Angle(1.,unit='deg'))
in_coords = coords.apply_region_mask(region)
out_coords = coords.apply_region_mask(region,inside=False)
assert in_coords.shape[0] == 1
assert in_coords.lon[0] == 0.
assert out_coords.lon[0] == 10.
@@ -6,6 +6,7 @@
from astropy.io import fits
from astropy.coordinates import SkyCoord
from astropy.units import Quantity
import astropy.units as u

This comment has been minimized.

@cdeil

cdeil May 29, 2018

Member

Remove unused import

@cdeil

cdeil May 29, 2018

Member

Remove unused import

This comment has been minimized.

@registerrier

registerrier Jun 1, 2018

Contributor

Done

@registerrier

registerrier Jun 1, 2018

Contributor

Done

from ..utils import fill_poisson
from ..geom import MapAxis, coordsys_to_frame
from ..base import Map
@@ -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,15 @@ 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 np.sum(maskmap.data) == 1

This comment has been minimized.

@cdeil

cdeil May 29, 2018

Member

Suggest to add an assert on the dtype of data. It's bool, no?

@cdeil

cdeil May 29, 2018

Member

Suggest to add an assert on the dtype of data. It's bool, no?

This comment has been minimized.

@registerrier

registerrier Jun 1, 2018

Contributor

Done

@registerrier

registerrier Jun 1, 2018

Contributor

Done

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,40 @@ 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
TODO: implement list of region for each axis

This comment has been minimized.

@cdeil

cdeil May 29, 2018

Member

I think in my last comment I suggested to remove this TODO (and also in the code below) and to instead show how to do this via a docs example.

@cdeil

cdeil May 29, 2018

Member

I think in my last comment I suggested to remove this TODO (and also in the code below) and to instead show how to do this via a docs example.

This comment has been minimized.

@registerrier

registerrier Jun 1, 2018

Contributor

Well it is not so easy to do actually. I will add an issue about this.

@registerrier

registerrier Jun 1, 2018

Contributor

Well it is not so easy to do actually. I will add an issue about this.

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
"""
from regions import PixCoord
# TODO : if Pixel Compound regions are taken into account, rather convert to PixelRegion

This comment has been minimized.

@cdeil

cdeil May 18, 2018

Member

Remove this TODO?
I don't understand it, and it doesn't seem to be needed?

@cdeil

cdeil May 18, 2018

Member

Remove this TODO?
I don't understand it, and it doesn't seem to be needed?

This comment has been minimized.

@registerrier

registerrier May 18, 2018

Contributor

Compound regions will be needed for exclusion masks, no?

@registerrier

registerrier May 18, 2018

Contributor

Compound regions will be needed for exclusion masks, no?

# if isinstance(region, SkyRegion):
# region = region.to_pixel(self.wcs)
if isinstance(region, SkyRegion):

This comment has been minimized.

@cdeil

cdeil May 29, 2018

Member

I think this would be simpler and equivalent to how regions works now:

if isinstance(region, SkyRegion):
    region = region.to_pixel(self.wcs)

pix = PixCoord(*self.get_idx())
return region.contains(pix)
@cdeil

cdeil May 29, 2018

Member

I think this would be simpler and equivalent to how regions works now:

if isinstance(region, SkyRegion):
    region = region.to_pixel(self.wcs)

pix = PixCoord(*self.get_idx())
return region.contains(pix)

This comment has been minimized.

@registerrier

registerrier Jun 1, 2018

Contributor

Done

@registerrier

registerrier Jun 1, 2018

Contributor

Done

coords=self.get_coord()
# Test to adapt to various call for region.contains
try:
res = region.contains(coords.skycoord, self.wcs)
except TypeError:
res = region.contains(coords.skycoord)
return res
else:
res = self.get_idx()
pcoords = PixCoord(res[0], res[1])
return region.contains(pcoords)
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)

This comment has been minimized.

@joleroi

joleroi May 25, 2018

Contributor

Can't you pass inside to geom.get_region_mask_array in order to prevent the duplicated check for if inside in the following lines?

@joleroi

joleroi May 25, 2018

Contributor

Can't you pass inside to geom.get_region_mask_array in order to prevent the duplicated check for if inside in the following lines?

This comment has been minimized.

@registerrier

registerrier May 25, 2018

Contributor

OK.

@registerrier
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)
ProTip! Use n and p to navigate between commits in a pull request.