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 WCS map cutout method #1446

Merged
merged 9 commits into from Jun 29, 2018
Merged

Add WCS map cutout method #1446

merged 9 commits into from Jun 29, 2018

Conversation

@AtreyeeS
Copy link
Member

@AtreyeeS AtreyeeS commented Jun 28, 2018

Hello,

I added a make_cutout() following the one in SkyCube.
This makes only a spatial cutout on all axis of the cube.

But in practise, I think the users will want to be able to make cutouts on all other axes as well, no?
I mean, maybe something like
cutout=cmap.make_cutout(position=pos,radius=rad, energy=(1,3)*u.TeV, time=(t1,t2) )

Or can cutout remain only for spatial cuts, and then other axes extracted as needed?

@cdeil cdeil added the feature label Jun 28, 2018
@cdeil cdeil added this to the 0.8 milestone Jun 28, 2018
Copy link
Member

@cdeil cdeil left a comment

@AtreyeeS - I left some quick inline comments.

Suggest to just do the spatial cutout in this PR, and discuss whether we want to combine slicing with other axes in the call tomorrow.

So this already works if you have extra non-spatial axes? If yes, please add a test case, or change the existing one to have one or two small extra axes.

@adonath - I'll leave the rest of the review and putting this in to you.

margin : `~astropy.coordinates.Angle`, optional
Additional safety margin. If specified must be same length as radius
mode : {'trim', 'partial', 'strict'}
Mode option for Cutout2D, see http://docs.astropy.org/en/stable/api/astropy.nddata.utils.Cutout2D.html

This comment has been minimized.

@cdeil

cdeil Jun 28, 2018
Member

Use this instead of a URL to link to the Astropy docs:

`~astropy.nddata.utils.Cutout2D`

This is using intersphinx to make the links.

"""

from astropy.nddata import Cutout2D

This comment has been minimized.

@cdeil

cdeil Jun 28, 2018
Member

Move import to the top of the file.


from astropy.nddata import Cutout2D

size=radius+margin if margin else radius

This comment has been minimized.

@cdeil

cdeil Jun 28, 2018
Member

Write PEP8 clean code, i.e. have a space on the left and right of =.
Use an editor (like PyCharm or VS Code with Python extension) that does that automatically or via "format code" for you; using command line tools like autopep8 is too painful.

position=position, size=size, mode=mode)

# Create the slices with the non-spatial axis
cs=()

This comment has been minimized.

@cdeil

cdeil Jun 28, 2018
Member

Use a list and apped isntead of a tuple.
Tuples are immutable and you keep making copies and new ones, where really you want to a mutable data structure you can append to here -> a list.

geom = WcsGeom(cutout2d.wcs, cutout2d.shape[::-1], axes=self.geom.axes)
data = self.data[cutout_slices]

ndmap_cutout = WcsNDMap(geom, data)

This comment has been minimized.

@cdeil

cdeil Jun 28, 2018
Member

This local variable isn't useful, just

return WcsNDMap(geom, data)
@cdeil cdeil added this to To do in gammapy.maps via automation Jun 28, 2018
@cdeil cdeil changed the title Added make_cutout() in wcsnd.py Add WCS map cutout method Jun 28, 2018
@AtreyeeS
Copy link
Member Author

@AtreyeeS AtreyeeS commented Jun 28, 2018

Yes, this works for any number of non-spatial axes. Changed the test to make a cutout on a 4d map now (and other inline corrections implemented)

position : `~astropy.coordinates.SkyCoord`
Center position of the cutout box
radius : tuple of `~astropy.coordinates.Angle`

This comment has been minimized.

@adonath

adonath Jun 28, 2018
Member

I think the radius argument is misleading, because rectangular cut-outs are possible. I would suggest to rename it to width and also document the axis order (lon, lat) for the user.

This comment has been minimized.

Center position of the cutout box
radius : tuple of `~astropy.coordinates.Angle`
Angular sizes of the box
margin : `~astropy.coordinates.Angle`, optional

This comment has been minimized.

@adonath

adonath Jun 28, 2018
Member

Can you shortly explain the motivation for the margin parameter? As the user could just compute radius = radius + margin outside the make_cutout function with the same result...so I would suggest to remove it.

This comment has been minimized.

@AtreyeeS

AtreyeeS Jun 28, 2018
Author Member

The margin keyword was there in the SkyCube implementation of cutout. I guess it was to ensure an overlap while making mosaics for example...
But I would be happy to remove it!

-------
cutout : `~gammapy.maps.WcsNDMap`
The cutout map itself
cutout_slices : tuple

This comment has been minimized.

@adonath

adonath Jun 28, 2018
Member

Remove the cutout_slices from the docstring, as it is not returned.

proj='CAR', coordsys='GAL')
m = WcsNDMap(geom, data=np.ones((10, 10)), unit='m2')
pos=SkyCoord(0, 0, unit='deg', frame='galactic')
cutout=m.make_cutout(position=pos,radius=3.0*u.deg,margin=0.1*u.deg)

This comment has been minimized.

@adonath

adonath Jun 28, 2018
Member

Can you add a test case with a rectangular cut-out e.g. (3 * u.deg, 2 * u.deg) and make sure the axis order and data shape of the returned map is correct?

cs=cs+tuple([slice(0,self.data.shape[i])])
cutout_slices=tuple(cs)+cutout2d.slices_original

geom = WcsGeom(cutout2d.wcs, cutout2d.shape[::-1], axes=self.geom.axes)

This comment has been minimized.

@adonath

adonath Jun 28, 2018
Member

I would suggest to move the part of the code that acts on the WcsGeom object to a WcsGeom.cutout() method. But this is more of a general question. @cdeil @registerrier Do we want to mirror the methods, that don't act on data in the API of the MapGeom objects? I guess it leads to simpler implementations but also duplicates some documentation...but let's discuss this tomorrow in the call. It should not stop this PR.

geom = WcsGeom(cutout2d.wcs, cutout2d.shape[::-1], axes=self.geom.axes)
data = self.data[cutout_slices]

ndmap_cutout = WcsNDMap(geom, data)

This comment has been minimized.

@adonath

adonath Jun 28, 2018
Member

Note that one should pass the other map attributes such as meta and unit as well.

cs=()
for i in range(self.data.ndim-2):
cs=cs+tuple([slice(0,self.data.shape[i])])
cutout_slices=tuple(cs)+cutout2d.slices_original

This comment has been minimized.

@adonath

adonath Jun 28, 2018
Member

Note this can be achieved much simpler using the Ellipsis object (that I actually learned about only few days ago).
Here is a minimal numpy example:

a = np.arange(1E6).reshape((10, 10, 100, 100))
print(a.shape)

idx = (Ellipsis, slice(25, 75), slice(25, 75))
b = a[idx]
print(b.shape)

This comment has been minimized.

@AtreyeeS

AtreyeeS Jun 28, 2018
Author Member

Thanks so much for this! It makes life so much simpler!

@adonath
Copy link
Member

@adonath adonath commented Jun 28, 2018

@AtreyeeS Thanks for this PR! Please note my review comments (some of them are marked as outdated by github, but might still be open to address)

AtreyeeS added 3 commits Jun 28, 2018
@AtreyeeS
Copy link
Member Author

@AtreyeeS AtreyeeS commented Jun 28, 2018

Thanks @adonath,

Made the changes as you suggested.
Have not made any additions to the WcsGeom class for now.. lets discuss this tomorrow.

Copy link
Member

@adonath adonath left a comment

I've made a few final cosmetic changes. Once Travis-CI has passed I'll merge the PR

@registerrier registerrier moved this from To do to In progress in gammapy.maps Jun 29, 2018
idx = (0,) * len(self.geom.axes)

cutout2d = Cutout2D(data=self.data[idx], wcs=self.geom.wcs,
position=position, size=width, mode=mode, copy=copy)

This comment has been minimized.

@adonath

adonath Jun 29, 2018
Member

Note that the data passed to Cutout2D is not actually used, so it does not help make the copy here. Instead just add

if copy:
    data = data.copy()

to the end of the cutout method.

@adonath adonath merged commit 5afddfc into gammapy:master Jun 29, 2018
1 of 2 checks passed
1 of 2 checks passed
continuous-integration/travis-ci/pr The Travis CI build failed
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
gammapy.maps automation moved this from In progress to Done Jun 29, 2018
@AtreyeeS AtreyeeS deleted the AtreyeeS:cutout branch Jul 12, 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