Skip to content
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
Merged

Add map region mask #1421

merged 7 commits into from Jun 3, 2018

Conversation

@registerrier
Copy link
Contributor

@registerrier registerrier commented May 18, 2018

This PR adds a make_region_maskmethod to HpxNDMap and WcsNDMap.

@registerrier registerrier added this to To do in gammapy.maps via automation May 18, 2018
@registerrier registerrier self-assigned this May 18, 2018
@registerrier registerrier added this to the 0.8 milestone May 18, 2018
@registerrier registerrier requested a review from cdeil May 18, 2018
Copy link
Member

@cdeil cdeil left a comment

@registerrier - Thanks!

I left some inline comments.

My main comment / question is concerning if / how to do contains for sky regions. I think the discussions in astropy-regions have settled on the understanding that sky regions contains isn't well-defined in general and not supported, one always needs a WCS and then does the contains in pixel coordinates, like this:
http://docs.gammapy.org/dev/_modules/gammapy/image/core.html#SkyImage.region_mask

The exception is circle and polygon regions, for those a spherical contains is well-defined, and those should be implemented in astropy-regions as special classes. But that's not there, so I don't think at the moment we should support it here, which makes we don't support region filtering for HPX at all.

@joleroi @registerrier @woodmd - Thoughts?

@@ -1731,6 +1732,23 @@ def solid_angle(self):
import healpy as hp
return Quantity(hp.nside2pixarea(self.nside), 'sr')

def get_region_idx(self, region):
"""Return idx of pixels inside region

This comment has been minimized.

@cdeil

cdeil May 18, 2018
Member

I don't think this description is correct. It returns a boolean mask, not a pixel idx array, no?

@@ -648,3 +649,28 @@ def get_angle(x, t):
ax.coords.grid(color='w', linestyle=':', linewidth=0.5)

return fig, ax, p

def make_region_mask(self, region, inside=True):

This comment has been minimized.

@cdeil

cdeil May 18, 2018
Member

I think this code is 100% identical in WCS and HPX case.
Suggest to avoid the code duplication and move it to the base class.

This comment has been minimized.

@registerrier

registerrier May 18, 2018
Author Contributor

A priori no because you have no constructor on the base class.
I suppose we can do it using Map.from_geom()and then assign Map.data=mask

if inside is False:
np.logical_not(mask,out=mask)

# TODO : update meta table to include something about the region used for mask creation?

This comment has been minimized.

@cdeil

cdeil May 18, 2018
Member

Suggest to remove this TODO.

Adding such book-keeping could be done all throughout Gammapy, but would be more than double of actually writing the actual code and needs a lot of thought. Let's write code first, and then think about such provenance later.

This comment has been minimized.

@registerrier

registerrier May 18, 2018
Author Contributor

OK will remove the TODO.
I put it to remember that at some point we need to do something here with the metadata.

@@ -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=(11,11), binsz=2,

This comment has been minimized.

@cdeil

cdeil May 18, 2018
Member

Suggest to do fewer pixels, so that printouts are nice if an error occurs.
This is what we had before:

def test_region_mask():

Another key question is if we return dtype of bool or int. Please add an assert to establish that.
I think bool makes sense, but then write fails because FITS can't handle it? If so, we could consider changing write to convert to unit8 for bool maps.

This comment has been minimized.

@cdeil

cdeil May 18, 2018
Member

Another thing you could do for the test is to put it near the lon = 0 / 360 line, to see what happens on the other side and establish behaviour. This was a cause of bugs in the HESS software for many years, so it would help to at least establish behaviour or to get it completely right from the start in Gammapy.

"""

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?

This comment has been minimized.

@registerrier

registerrier May 18, 2018
Author Contributor

Compound regions will be needed for exclusion masks, no?

@cdeil cdeil changed the title Map mask Add map region mask May 18, 2018
def make_region_mask(self, region, inside=True):
"""Create a mask of a given region
TODO: implement list of region for each axis

This comment has been minimized.

@cdeil

cdeil May 18, 2018
Member

Concerning this TODO:

How should we handle extra axes?

I'm -1 on trying to define ND regions, that's a huge project.

That leaves the option of making a 2D mask, or the same 2D mask duplicated along the axes.

Here we did the replicated along the energy axis:
http://docs.gammapy.org/dev/_modules/gammapy/cube/core.html#SkyCube.region_mask
I think for solid_angle that's also what you did @registerrier ?
I would be in favour of returning 2D by default, but adding a bool option which replicates and makes an ND mask.

And then I'd also suggest to add an example how to do something more fancy, e.g. if there is an energy-dependent FOV one can make an empty mask, then call iter_by_image and fill as you like:
http://docs.gammapy.org/dev/maps/index.html#iterating-on-a-map
Inventing an API for that is difficult, and things can be added later, so for now I'd suggest to solve this via docs.

This comment has been minimized.

@registerrier

registerrier May 18, 2018
Author Contributor

Well that was the idea basically.
You can provide a list/tuple of SkyRegions and then apply iteration.

@joleroi
Copy link
Contributor

@joleroi joleroi commented May 18, 2018

I don't think at the moment we should support it here, which makes we don't support region filtering for HPX at all.

Agreed, let's leave this to astropy.regions

raise NotImplementedError

coords=self.get_coord()
return region.contains(coords.skycoord)

This comment has been minimized.

@joleroi

joleroi May 18, 2018
Contributor

Note that in future regions releases contains will require a wcs object (see latest docs). This is due to the fact, that SkyRegions are basically PixelRegions with the coordinates defined on the sky (as opposed to real sky regions).

So this method should work for PixelRegions - at the moment you raise a NotImplementedError on that. For convenience you could accept SkyRegion plus a wcs object. But like it is now this method will stop working soon.

----------
region : `~regions.SkyRegion` object
A region on the sky could be defined in pixel or sky coordinates.
Note that `~regions.PixelRegion` does not work for healpix maps.

This comment has been minimized.

@joleroi

joleroi May 18, 2018
Contributor

As mentioned by @cdeil we cannot really support HPX at the moment

This comment has been minimized.

@registerrier

registerrier May 18, 2018
Author Contributor

OK will remove it then.

@registerrier
Copy link
Contributor Author

@registerrier registerrier commented May 24, 2018

The new commits have:

  • removed the region filtering in map.HpxGeom and its test function
  • changed the name of the method toWcsGeom.get_region_mask_array
  • include the new calling scheme for SkyRegion.contains
  • added a MapCoord.apply_region_mask and its associated test

Merge?

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?

This comment has been minimized.

@registerrier

registerrier May 25, 2018
Author Contributor

OK.

@registerrier
Copy link
Contributor Author

@registerrier registerrier commented May 25, 2018

I realize that there is a serious issue with MapCoord.apply_region_mask. At the moment is uses the simple SkyRegion.contains(skycoord) method which is about to be replaced by SkyRegion.contains(skycoord, wcs).
There is no WCS associated to a MapCoord, so we have to pass one to make the selection.

This is a serious issue with regionsno? A SkyRegion should be a region on the sphere no? Otherwise you end up with different selections depending on the WCS object you pass...

@joleroi
Copy link
Contributor

@joleroi joleroi commented May 25, 2018

This is a serious issue with regionsno? A SkyRegion should be a region on the sphere no? Otherwise you end up with different selections depending on the WCS object you pass...

Yes, to my understanding this is intended for SkyRegion. It is really a region on some kind of sky projection. The plan is to add SphericalSkyRegion later, to cover the use case were people want to put a region on an actual sphere. SphericalSkyRegion will implement circles and polygons. There is some discussion on this in astropy/regions#107

Copy link
Member

@cdeil cdeil left a comment

@registerrier - I left some more small inline comments.

Do you have time to wrap this up?
(feel free to merge, no more iteration needed with me)

@@ -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

This comment has been minimized.

@registerrier

registerrier Jun 1, 2018
Author Contributor

Done

@@ -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

This comment has been minimized.

@registerrier

registerrier Jun 1, 2018
Author Contributor

Done

@@ -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

This comment has been minimized.

@registerrier

registerrier Jun 1, 2018
Author Contributor

Done

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.

This comment has been minimized.

@registerrier

registerrier Jun 1, 2018
Author Contributor

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

# 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)

This comment has been minimized.

@registerrier

registerrier Jun 1, 2018
Author Contributor

Done

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?

This comment has been minimized.

@registerrier

registerrier Jun 1, 2018
Author Contributor

Done

@registerrier registerrier moved this from To do to In progress in gammapy.maps Jun 1, 2018
@registerrier
Copy link
Contributor Author

@registerrier registerrier commented Jun 1, 2018

OK I have made the modifications.
There are some difficulties which I had not foreseen with multi resolution maps. I have added an Exception if one try to use one. This will need to be sorted out later.
I have also removed for the moment the MapCoord region filtering. This might be better handled with a MapMask for instance.

@registerrier registerrier merged commit 8ab962f into gammapy:master Jun 3, 2018
1 of 2 checks passed
1 of 2 checks passed
continuous-integration/travis-ci/pr The Travis CI build could not complete due to an error
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
gammapy.maps automation moved this from In progress to Done Jun 3, 2018
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
gammapy.maps
  
Done
Linked issues

Successfully merging this pull request may close these issues.

None yet

3 participants