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

Remove old HEALPix image and cube classes #1169

Merged
merged 9 commits into from Oct 29, 2017

Conversation

Projects
None yet
3 participants
@cdeil
Member

cdeil commented Oct 17, 2017

@adonath - I'd like to remove the SkyImageHealpix and SkyCubeHealpix classes that you added in January (see #841).

My understanding is that these classes were never exposed via the Gammapy docs, except for the one caller of SkyCubeHealpix at http://docs.gammapy.org/en/latest/api/gammapy.datasets.FermiLATDataset.html#gammapy.datasets.FermiLATDataset.exposure and that currently the only use case is exactly that case of reading that HEALPix exposure cube and reprojecting it onto a WCS cube?

I think if that is correct, then it should be easy to switch to HpxMapND and keep that use case fully working (and the sooner we remove it, the less likely is that others will start using it and their script breaks later on).

@adonath - If I keep TestFermiLATDataset at

class TestFermiLATDataset:
as well as https://nbviewer.jupyter.org/github/gammapy/gammapy-extra/blob/master/notebooks/data_fermi_lat.ipynb#Exposure working, are you OK with this change? Or are there other use-cases I should watch out for? If yes, can you please leave a code snippet with a minimal example that I could add as a test case?

@cdeil cdeil added the cleanup label Oct 17, 2017

@cdeil cdeil added this to the 0.7 milestone Oct 17, 2017

@cdeil cdeil requested a review from adonath Oct 17, 2017

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Oct 17, 2017

Member

@woodmd - I'm trying to read the HEALPix counts and exposure cube with gammapy.maps, that was written by a (probably a bit outdated) version of the Fermi science tools. See the files fermi_2fhl_counts_cube_hpx.fits.gz and fermi_2fhl_exposure_cube_hpx.fits.gz here:
https://github.com/gammapy/gammapy-fermi-lat-data/tree/master/2fhl

I'm getting this error:

>>> from gammapy.maps import HpxMapND
>>> HpxMapND.read('https://github.com/gammapy/gammapy-fermi-lat-data/raw/master/2fhl/fermi_2fhl_counts_cube_hpx.fits.gz')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/deil/code/gammapy/gammapy/maps/base.py", line 126, in read
    map_out = cls.from_hdulist(hdulist, **kwargs)
  File "/Users/deil/code/gammapy/gammapy/maps/hpxmap.py", line 116, in from_hdulist
    return cls.from_hdu(hdu, hdu_bands)
  File "/Users/deil/code/gammapy/gammapy/maps/hpxcube.py", line 74, in from_hdu
    if hdu.header['INDXSCHM'] == 'SPARSE':
  File "/Users/deil/Library/Python/3.6/lib/python/site-packages/astropy/io/fits/header.py", line 145, in __getitem__
    card = self._cards[self._cardindex(key)]
  File "/Users/deil/Library/Python/3.6/lib/python/site-packages/astropy/io/fits/header.py", line 1651, in _cardindex
    raise KeyError("Keyword {!r} not found.".format(keyword))
KeyError: "Keyword 'INDXSCHM' not found."

Here you explain about the conv option for write:
http://docs.gammapy.org/en/latest/maps/index.html#fits-i-o
Maybe the same or a similar option is needed for read?
Or the read method can be made smarter to just work?
Or I should read it in some other way that doesn't require this header key?

Member

cdeil commented Oct 17, 2017

@woodmd - I'm trying to read the HEALPix counts and exposure cube with gammapy.maps, that was written by a (probably a bit outdated) version of the Fermi science tools. See the files fermi_2fhl_counts_cube_hpx.fits.gz and fermi_2fhl_exposure_cube_hpx.fits.gz here:
https://github.com/gammapy/gammapy-fermi-lat-data/tree/master/2fhl

I'm getting this error:

>>> from gammapy.maps import HpxMapND
>>> HpxMapND.read('https://github.com/gammapy/gammapy-fermi-lat-data/raw/master/2fhl/fermi_2fhl_counts_cube_hpx.fits.gz')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/deil/code/gammapy/gammapy/maps/base.py", line 126, in read
    map_out = cls.from_hdulist(hdulist, **kwargs)
  File "/Users/deil/code/gammapy/gammapy/maps/hpxmap.py", line 116, in from_hdulist
    return cls.from_hdu(hdu, hdu_bands)
  File "/Users/deil/code/gammapy/gammapy/maps/hpxcube.py", line 74, in from_hdu
    if hdu.header['INDXSCHM'] == 'SPARSE':
  File "/Users/deil/Library/Python/3.6/lib/python/site-packages/astropy/io/fits/header.py", line 145, in __getitem__
    card = self._cards[self._cardindex(key)]
  File "/Users/deil/Library/Python/3.6/lib/python/site-packages/astropy/io/fits/header.py", line 1651, in _cardindex
    raise KeyError("Keyword {!r} not found.".format(keyword))
KeyError: "Keyword 'INDXSCHM' not found."

Here you explain about the conv option for write:
http://docs.gammapy.org/en/latest/maps/index.html#fits-i-o
Maybe the same or a similar option is needed for read?
Or the read method can be made smarter to just work?
Or I should read it in some other way that doesn't require this header key?

@adonath

This comment has been minimized.

Show comment
Hide comment
@adonath

adonath Oct 17, 2017

Member

@cdeil This is exactly the only use case I had for the SkyCubeHealpix class; reading the Fermi-LAT exposure cube and reproject it to WCS. As long as this works with the new maps classes I'm fine to remove those old classes.

Note that this might also require a change in the FermiLATBasicImageEstimator class, because the reprojected cube must be a SkyCube again. This should be covered with the tests in gammapy/image/tests/test_basic.py but only locally, because it requires the $GAMMAPY_FERMI_LAT_DATA to be available.

Member

adonath commented Oct 17, 2017

@cdeil This is exactly the only use case I had for the SkyCubeHealpix class; reading the Fermi-LAT exposure cube and reproject it to WCS. As long as this works with the new maps classes I'm fine to remove those old classes.

Note that this might also require a change in the FermiLATBasicImageEstimator class, because the reprojected cube must be a SkyCube again. This should be covered with the tests in gammapy/image/tests/test_basic.py but only locally, because it requires the $GAMMAPY_FERMI_LAT_DATA to be available.

@woodmd

This comment has been minimized.

Show comment
Hide comment
@woodmd

woodmd Oct 17, 2017

Member

The intention was that the read method would automatically detect the format convention of the file but it appears this isn't working. I'll have a closer look at this today.

Member

woodmd commented Oct 17, 2017

The intention was that the read method would automatically detect the format convention of the file but it appears this isn't working. I'll have a closer look at this today.

@woodmd

This comment has been minimized.

Show comment
Hide comment
@woodmd

woodmd Oct 17, 2017

Member

I made some small changes to support reading files that are missing INDXSCHM and/or COORDSYS header keywords. The current logic for reading non-standard formats is a bit of a mess so we should think about adding some explicit tests in gammapy.maps for reading files generated with the Fermi STs.

Member

woodmd commented Oct 17, 2017

I made some small changes to support reading files that are missing INDXSCHM and/or COORDSYS header keywords. The current logic for reading non-standard formats is a bit of a mess so we should think about adding some explicit tests in gammapy.maps for reading files generated with the Fermi STs.

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Oct 17, 2017

Member

@woodmd - Thanks for fixing the read!

It's working now. There's still this code that needs to be adapted for doing the reproject:
https://gist.github.com/cdeil/19762df1038889985744071c376f9281

Another question that comes up is unit handling. Before we had a Quantity object for data. Now it's a numpy array and the unit info isn't read. I think I'll leave like that for now, and then add unit to the gammapy.maps classes in a separate PR later.

Member

cdeil commented Oct 17, 2017

@woodmd - Thanks for fixing the read!

It's working now. There's still this code that needs to be adapted for doing the reproject:
https://gist.github.com/cdeil/19762df1038889985744071c376f9281

Another question that comes up is unit handling. Before we had a Quantity object for data. Now it's a numpy array and the unit info isn't read. I think I'll leave like that for now, and then add unit to the gammapy.maps classes in a separate PR later.

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Oct 18, 2017

Member

I've done some docs changes in 2e86dec , mainly to SkyImage and SkyCube:

  • start the pattern of linking from class-level docstrings to the high-level docs page
  • start the pattern of linking from high-level docs page to tutorials. For now I've put it in the "using" section, instead of at the end of the "getting started" section like @joleroi did in gammapy.spectrum before, but I don't feel strongly either way. (my reasoning was that like this it's easier to share a link to more info on "using").
  • add a note encouraging people to move from SkyImage and SkyCube to use gammapy.maps. I'm not making it a warning / deprecation yet: this needs further discussion. But pointing them to gammapy.maps and letting them know that this is coming is useful I think.

To make the "reproject" work here, I will now add converter helper methods to SkyImage and SkyCube, to convert to the corresponding gammapy.maps map and geometry objects. I'm not sure yet if I should go "all in" with gammapy.maps for the Fermi dataset and basic image estimator. At the moment it's a mix, we're still using SkyImage and SkyCube, but for HEALPix the gammapy.maps class.

@adonath - Thoughts?

Member

cdeil commented Oct 18, 2017

I've done some docs changes in 2e86dec , mainly to SkyImage and SkyCube:

  • start the pattern of linking from class-level docstrings to the high-level docs page
  • start the pattern of linking from high-level docs page to tutorials. For now I've put it in the "using" section, instead of at the end of the "getting started" section like @joleroi did in gammapy.spectrum before, but I don't feel strongly either way. (my reasoning was that like this it's easier to share a link to more info on "using").
  • add a note encouraging people to move from SkyImage and SkyCube to use gammapy.maps. I'm not making it a warning / deprecation yet: this needs further discussion. But pointing them to gammapy.maps and letting them know that this is coming is useful I think.

To make the "reproject" work here, I will now add converter helper methods to SkyImage and SkyCube, to convert to the corresponding gammapy.maps map and geometry objects. I'm not sure yet if I should go "all in" with gammapy.maps for the Fermi dataset and basic image estimator. At the moment it's a mix, we're still using SkyImage and SkyCube, but for HEALPix the gammapy.maps class.

@adonath - Thoughts?

@adonath

This comment has been minimized.

Show comment
Hide comment
@adonath

adonath Oct 18, 2017

Member

I agree it is a good idea to add the converter methods to the existing SkyImage and SkyCube methods; this will help a lot with the transition to the new map classes. But in general think it is still too early to go "all in" with the new map classes, even for the FermiLATBasicImageEstimator. E.g. some important methods that are missing right now are .cutout(), .paste() and .convolve(), which are not completely trivial to implement. From my side it is completely fine for now to use the new map classes for healpix and keep the old ones for WCS based images.

Member

adonath commented Oct 18, 2017

I agree it is a good idea to add the converter methods to the existing SkyImage and SkyCube methods; this will help a lot with the transition to the new map classes. But in general think it is still too early to go "all in" with the new map classes, even for the FermiLATBasicImageEstimator. E.g. some important methods that are missing right now are .cutout(), .paste() and .convolve(), which are not completely trivial to implement. From my side it is completely fine for now to use the new map classes for healpix and keep the old ones for WCS based images.

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Oct 18, 2017

Member

I've added SkyCube.from_wcs_map_nd and SkyCube.to_wcs_map_nd helper functions in 8dc98a5 and a minimal test case.

@woodmd - One mistake I made at first was to pass an astropy.wcs.WCS object to gammapy.maps.WcsMapND, because the first parameter is called wcs. I think maybe renaming to geom or even wcs_geom could be helpful, i.e. only use "wcs" for WCS objects?

Another thing that wasn't obvious to me from http://docs.gammapy.org/en/latest/maps/index.html was the axis order: in SkyCube and 3D WCS-based FITS files we have the data array with axis order energy, lat, lon. In the text you say "By convention the first two elements in the coordinate tuple correspond to longitude and latitude followed by one element for each non-spatial dimension." but this is pretty hard to find and it's still not clear if energy can come before the spatial axes for the data I pass in to WcsMapNd or not. I'm not quite sure what docs additions would help there. Maybe just showing a "print-out" that shows the axis order of the data array like at http://docs.gammapy.org/en/latest/cube/index.html#getting-started would be enough, maybe a separate section to discuss axis order would be useful?

So now the TestFermiLATBasicImageEstimator.test_run test fails like this:
https://gist.github.com/cdeil/cf5c5d39dd87e7c4fe0633bf378d0e61
because your reproject doesn't automatically also do interpolation for the energy axis.
I think I have to do a get_by_coords call first: http://docs.gammapy.org/en/latest/maps/index.html#interpolation
I'll try that now. But @woodmd - again, from the high-level docs it's not obvious to me what I should pass if I just want to interpolate to a new energy grid, and keep the spatial grid unchanged. If you have time - I think adding a high-level docs example explaining how to do that would be helpful.

Member

cdeil commented Oct 18, 2017

I've added SkyCube.from_wcs_map_nd and SkyCube.to_wcs_map_nd helper functions in 8dc98a5 and a minimal test case.

@woodmd - One mistake I made at first was to pass an astropy.wcs.WCS object to gammapy.maps.WcsMapND, because the first parameter is called wcs. I think maybe renaming to geom or even wcs_geom could be helpful, i.e. only use "wcs" for WCS objects?

Another thing that wasn't obvious to me from http://docs.gammapy.org/en/latest/maps/index.html was the axis order: in SkyCube and 3D WCS-based FITS files we have the data array with axis order energy, lat, lon. In the text you say "By convention the first two elements in the coordinate tuple correspond to longitude and latitude followed by one element for each non-spatial dimension." but this is pretty hard to find and it's still not clear if energy can come before the spatial axes for the data I pass in to WcsMapNd or not. I'm not quite sure what docs additions would help there. Maybe just showing a "print-out" that shows the axis order of the data array like at http://docs.gammapy.org/en/latest/cube/index.html#getting-started would be enough, maybe a separate section to discuss axis order would be useful?

So now the TestFermiLATBasicImageEstimator.test_run test fails like this:
https://gist.github.com/cdeil/cf5c5d39dd87e7c4fe0633bf378d0e61
because your reproject doesn't automatically also do interpolation for the energy axis.
I think I have to do a get_by_coords call first: http://docs.gammapy.org/en/latest/maps/index.html#interpolation
I'll try that now. But @woodmd - again, from the high-level docs it's not obvious to me what I should pass if I just want to interpolate to a new energy grid, and keep the spatial grid unchanged. If you have time - I think adding a high-level docs example explaining how to do that would be helpful.

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Oct 18, 2017

Member

@adonath - Would it be OK if I change FermiLATBasicImageEstimator to take the energy array to use for the computation as input parameter?

At the moment the user has no control over it, it's taken from the diffuse background cube.

I think letting the user choose the energy grid for the computation, and interpolating the exposure and diffuse background (or in the future other model components) to that energy grid is simpler / more flexible, no?

Member

cdeil commented Oct 18, 2017

@adonath - Would it be OK if I change FermiLATBasicImageEstimator to take the energy array to use for the computation as input parameter?

At the moment the user has no control over it, it's taken from the diffuse background cube.

I think letting the user choose the energy grid for the computation, and interpolating the exposure and diffuse background (or in the future other model components) to that energy grid is simpler / more flexible, no?

@adonath

This comment has been minimized.

Show comment
Hide comment
@adonath

adonath Oct 18, 2017

Member

@cdeil Internally the FermiLATBasicImageEstimator uses the compute_npred_cube utility function, where the parameter for precision of the energy integration is integral_resolution right now. I think I would be nice to have this consistent.

I remember that I've thought about the question on which energy and spatial grid to compute the predicted number of counts from the background model. I agree it would be simplest to just reproject to the spatial and energy binning given by the user, but because of performance (and memory!) issues I decided to use the spatial and energy binning of the diffuse background model and only reproject in the latest step. The main reason is that it uses the available information in the background model in an optimal way. With any other choice of a different spatial or energy binning you will either waste computing time or throw away information.
While the first is certainly not desired, that latter might be interesting for users in some cases, where performance is really important. So I'm not completely against introducing the energy binning as an additional option, but I think the use cases for this are very limited and the current behaviour should be kept as the default.

Member

adonath commented Oct 18, 2017

@cdeil Internally the FermiLATBasicImageEstimator uses the compute_npred_cube utility function, where the parameter for precision of the energy integration is integral_resolution right now. I think I would be nice to have this consistent.

I remember that I've thought about the question on which energy and spatial grid to compute the predicted number of counts from the background model. I agree it would be simplest to just reproject to the spatial and energy binning given by the user, but because of performance (and memory!) issues I decided to use the spatial and energy binning of the diffuse background model and only reproject in the latest step. The main reason is that it uses the available information in the background model in an optimal way. With any other choice of a different spatial or energy binning you will either waste computing time or throw away information.
While the first is certainly not desired, that latter might be interesting for users in some cases, where performance is really important. So I'm not completely against introducing the energy binning as an additional option, but I think the use cases for this are very limited and the current behaviour should be kept as the default.

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Oct 18, 2017

Member

@adonath - OK, I'll do some limited changes here for now to make it work.

But just to note for the future: these image estimator classes really are cube followed by energy integration. And we just need cube analysis, where the diffuse background model is just one model component. Sometimes it has very coarse, sometimes very fine energy binning. It might not be appropriate for the user's analysis, e.g. a line search at 100 GeV.

So really, I think we should just allow the user to choose the energy and spatial binning, and to interpolate / integrate IRFs and background models onto that grid. Adopting energy or spatial grids from inputs like IRFs and background models isn't a good solution for the general case, and we don't really want to maintain multiple analysis codes that basically do the same thing.

But OK, I don't have time to write a better cube analysis and then on top of that image maker at the moment. So I'll try to make minimal changes here to get this existing Fermi image maker to run again and to give the same output it did before.

Member

cdeil commented Oct 18, 2017

@adonath - OK, I'll do some limited changes here for now to make it work.

But just to note for the future: these image estimator classes really are cube followed by energy integration. And we just need cube analysis, where the diffuse background model is just one model component. Sometimes it has very coarse, sometimes very fine energy binning. It might not be appropriate for the user's analysis, e.g. a line search at 100 GeV.

So really, I think we should just allow the user to choose the energy and spatial binning, and to interpolate / integrate IRFs and background models onto that grid. Adopting energy or spatial grids from inputs like IRFs and background models isn't a good solution for the general case, and we don't really want to maintain multiple analysis codes that basically do the same thing.

But OK, I don't have time to write a better cube analysis and then on top of that image maker at the moment. So I'll try to make minimal changes here to get this existing Fermi image maker to run again and to give the same output it did before.

@woodmd

This comment has been minimized.

Show comment
Hide comment
@woodmd

woodmd Oct 20, 2017

Member

One mistake I made at first was to pass an astropy.wcs.WCS object to gammapy.maps.WcsMapND, because the first parameter is called wcs. I think maybe renaming to geom or even wcs_geom could be helpful, i.e. only use "wcs" for WCS objects?

Fine with me. I can combine this with your other proposed change to add factory methods to instantiate an empty map.

Another thing that wasn't obvious to me from http://docs.gammapy.org/en/latest/maps/index.html was the axis order

Longitude and latitude are required to be the first and second dimension. However the data array and accessor methods use opposite conventions for axis ordering. Following the convention of FITS images the data array has a column-wise ordering such that the spatial dimensions are last. However accessor methods follow the row-wise ordering convention as used in astropy.wcs whereby longitude and latitude come first followed by non-spatial dimensions. We might think about whether it would be simpler to use row-wise ordering everywhere. For now I can try to improve the documentation to explain this a bit better.

So now the TestFermiLATBasicImageEstimator.test_run test fails like this:

Interpolation on non-spatial dimensions in not yet implemented for healpix maps. If this is the only thing blocking this PR I can move this up my priority list. With respect to the syntax, you would simply create a map geometry object with the desired binning in energy and pass that to the reproject method.

Member

woodmd commented Oct 20, 2017

One mistake I made at first was to pass an astropy.wcs.WCS object to gammapy.maps.WcsMapND, because the first parameter is called wcs. I think maybe renaming to geom or even wcs_geom could be helpful, i.e. only use "wcs" for WCS objects?

Fine with me. I can combine this with your other proposed change to add factory methods to instantiate an empty map.

Another thing that wasn't obvious to me from http://docs.gammapy.org/en/latest/maps/index.html was the axis order

Longitude and latitude are required to be the first and second dimension. However the data array and accessor methods use opposite conventions for axis ordering. Following the convention of FITS images the data array has a column-wise ordering such that the spatial dimensions are last. However accessor methods follow the row-wise ordering convention as used in astropy.wcs whereby longitude and latitude come first followed by non-spatial dimensions. We might think about whether it would be simpler to use row-wise ordering everywhere. For now I can try to improve the documentation to explain this a bit better.

So now the TestFermiLATBasicImageEstimator.test_run test fails like this:

Interpolation on non-spatial dimensions in not yet implemented for healpix maps. If this is the only thing blocking this PR I can move this up my priority list. With respect to the syntax, you would simply create a map geometry object with the desired binning in energy and pass that to the reproject method.

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Oct 20, 2017

Member

Interpolation on non-spatial dimensions in not yet implemented for healpix maps. If this is the only thing blocking this PR I can move this up my priority list.

To finish up this PR I need to replace this:
https://github.com/gammapy/gammapy/pull/1169/files#diff-9ad09e56cce0d72a5c68615ab259e880L93

It's not actually doing interpolation on the energy axis, it's looping over energy bins and doing a 2-dim image HPX to WCS reprojection, to take the exposure cube from HPX to WCS and to then be able to use it for analysis (at the moment all our analysis is WCS based).

I'm not sure how to best finish this up. I don't think I want to start coding on gammapy.maps here, I'm not yet familiar enough with it that this is likely to yield something good. Maybe I just do the equivalent of this loop over energy bins for now here in the FermiLatBasicImageEstimator that Axel wrote, and then maybe this gets replaced later by gammapy.maps method calls.

I think you don't have to change priorities. Maybe documenting the axis order or opening a separate issue to discuss the possible change you mentioned should be high priority, since that will affect almost anyone that tries to use the maps.

Member

cdeil commented Oct 20, 2017

Interpolation on non-spatial dimensions in not yet implemented for healpix maps. If this is the only thing blocking this PR I can move this up my priority list.

To finish up this PR I need to replace this:
https://github.com/gammapy/gammapy/pull/1169/files#diff-9ad09e56cce0d72a5c68615ab259e880L93

It's not actually doing interpolation on the energy axis, it's looping over energy bins and doing a 2-dim image HPX to WCS reprojection, to take the exposure cube from HPX to WCS and to then be able to use it for analysis (at the moment all our analysis is WCS based).

I'm not sure how to best finish this up. I don't think I want to start coding on gammapy.maps here, I'm not yet familiar enough with it that this is likely to yield something good. Maybe I just do the equivalent of this loop over energy bins for now here in the FermiLatBasicImageEstimator that Axel wrote, and then maybe this gets replaced later by gammapy.maps method calls.

I think you don't have to change priorities. Maybe documenting the axis order or opening a separate issue to discuss the possible change you mentioned should be high priority, since that will affect almost anyone that tries to use the maps.

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Oct 20, 2017

Member

@woodmd - Looking at this:
http://docs.gammapy.org/en/latest/maps/index.html#classes
http://docs.gammapy.org/en/latest/_modules/gammapy/maps/hpxsparse.html#HpxMapSparse
http://docs.gammapy.org/en/latest/api/gammapy.maps.HpxMapSparse.html

I see some things that could maybe be improved:

Does HpxMapSparse have a "geom" attribute? In addition to a "hpx" attribute? What is the difference? Or is this a mistake in the class hierarchy or just in the way Sphinx shows attributes?

As mentioned earlier, I find the __init__ sub-optimal and would suggest to consider splitting out some things into factory methods. That's of course a matter of taste, I know that several other developers prefer to offer options in __init__ to construct objects in multiple ways (SkyCoord is an extreme example for this).

Concretely: there's no way to pass in a SparseArray object, and from the docs it's not visible at all that the data attribute is a SparseArray. So I would suggest to support passing in a SparseArray, and to also add SparseArray in __all__ so that it shows up in the docs. Of course there's always the tradeoff between overloading the docs with classes that many users will not need to know about, but in this case, I think some users, and Gammapy developers will have to learn about SparseArray and it's good to have it in the API docs.

I'm not sure where to best put the "make zero-filled data array" code. We had a corresponding discussion for SkyImage and SkyCube here: #513
Now in gammapy.maps there are "geometry" objects that know about the data array shape, so it becomes more attractive to put it in __init__, and the opportunity to put it as a factory method on the geometry object opens up. I think this would be OK (and make for a super-simple __init__:

geom = "geometry object"
data = geom.make_data(fill=0, dtype=np.int32)
map = Map(geom, data)

But also a factory method on the maps, or even keeping it in __init__ would be OK for me. @woodmd - If you decide to keep it in __init__: an extra sentence and an example in the class docstring explaining why there is a dtype option in addition to data would be helpful.

Member

cdeil commented Oct 20, 2017

@woodmd - Looking at this:
http://docs.gammapy.org/en/latest/maps/index.html#classes
http://docs.gammapy.org/en/latest/_modules/gammapy/maps/hpxsparse.html#HpxMapSparse
http://docs.gammapy.org/en/latest/api/gammapy.maps.HpxMapSparse.html

I see some things that could maybe be improved:

Does HpxMapSparse have a "geom" attribute? In addition to a "hpx" attribute? What is the difference? Or is this a mistake in the class hierarchy or just in the way Sphinx shows attributes?

As mentioned earlier, I find the __init__ sub-optimal and would suggest to consider splitting out some things into factory methods. That's of course a matter of taste, I know that several other developers prefer to offer options in __init__ to construct objects in multiple ways (SkyCoord is an extreme example for this).

Concretely: there's no way to pass in a SparseArray object, and from the docs it's not visible at all that the data attribute is a SparseArray. So I would suggest to support passing in a SparseArray, and to also add SparseArray in __all__ so that it shows up in the docs. Of course there's always the tradeoff between overloading the docs with classes that many users will not need to know about, but in this case, I think some users, and Gammapy developers will have to learn about SparseArray and it's good to have it in the API docs.

I'm not sure where to best put the "make zero-filled data array" code. We had a corresponding discussion for SkyImage and SkyCube here: #513
Now in gammapy.maps there are "geometry" objects that know about the data array shape, so it becomes more attractive to put it in __init__, and the opportunity to put it as a factory method on the geometry object opens up. I think this would be OK (and make for a super-simple __init__:

geom = "geometry object"
data = geom.make_data(fill=0, dtype=np.int32)
map = Map(geom, data)

But also a factory method on the maps, or even keeping it in __init__ would be OK for me. @woodmd - If you decide to keep it in __init__: an extra sentence and an example in the class docstring explaining why there is a dtype option in addition to data would be helpful.

@woodmd

This comment has been minimized.

Show comment
Hide comment
@woodmd

woodmd Oct 20, 2017

Member

It's not actually doing interpolation on the energy axis, it's looping over energy bins and doing a 2-dim image HPX to WCS reprojection, to take the exposure cube from HPX to WCS and to then be able to use it for analysis (at the moment all our analysis is WCS based).

Ah I see. At present there's not a clean way of doing this. One option would be to just manually interpolate at the given pixels, e.g.

geom = WcsGeom(npix=100)
m_hpx = HpxMapND.read('expsoure.fits')
for i, egy in enumerate(energies):
    m_wcs = WcsMapND(geom)
    coords = geom.get_coords()
    coords_hpx = tuple(list(coords) + [egy])
    vals = m_hpx.get_by_coords(coords_hpx, interp='linear')
    m_wcs.set_by_coord(coords, vals)

Obviously this is pretty cumbersome. The current reproject method always returns an object with the same dimensionality. To handle what you're doing here we could add an explicit method that reprojects a map to a map of lower dimensionality (e.g. cube to image), e.g. a new method called reproject_slice or reproject_to_slice. Alternatively we could add some new keywords to reproject that instruct the method to drop dimensions from the output map that are not defined in the input geometry.

Does HpxMapSparse have a "geom" attribute? In addition to a "hpx" attribute? What is the difference? Or is this a mistake in the class hierarchy or just in the way Sphinx shows attributes?

Yes it does but it points to the same object so we should probably just get rid of the hpx attribute.

So I would suggest to support passing in a SparseArray, and to also add SparseArray in all so that it shows up in the docs.

I think the current implementation of __init__ should already support passing in a SparseArray but the docstring should probably be updated to reflect that. I've also expanded the SparseAarray class documentation a bit so it probably makes sense to have it in __all__ now.

But also a factory method on the maps, or even keeping it in init would be OK for me.

I'd prefer having a single method to create an empty map object. Creating an empty map is a very common operation and requiring two function calls seems unnecessarily complex. In practice MapBase.create is a factory method that already partially serves this purpose -- however it doesn't accept a geometry argument but creates a geometry using the method arguments. Maybe we could move the logic in __init__ to a from_geom method that takes a geometry object and a data type.

Member

woodmd commented Oct 20, 2017

It's not actually doing interpolation on the energy axis, it's looping over energy bins and doing a 2-dim image HPX to WCS reprojection, to take the exposure cube from HPX to WCS and to then be able to use it for analysis (at the moment all our analysis is WCS based).

Ah I see. At present there's not a clean way of doing this. One option would be to just manually interpolate at the given pixels, e.g.

geom = WcsGeom(npix=100)
m_hpx = HpxMapND.read('expsoure.fits')
for i, egy in enumerate(energies):
    m_wcs = WcsMapND(geom)
    coords = geom.get_coords()
    coords_hpx = tuple(list(coords) + [egy])
    vals = m_hpx.get_by_coords(coords_hpx, interp='linear')
    m_wcs.set_by_coord(coords, vals)

Obviously this is pretty cumbersome. The current reproject method always returns an object with the same dimensionality. To handle what you're doing here we could add an explicit method that reprojects a map to a map of lower dimensionality (e.g. cube to image), e.g. a new method called reproject_slice or reproject_to_slice. Alternatively we could add some new keywords to reproject that instruct the method to drop dimensions from the output map that are not defined in the input geometry.

Does HpxMapSparse have a "geom" attribute? In addition to a "hpx" attribute? What is the difference? Or is this a mistake in the class hierarchy or just in the way Sphinx shows attributes?

Yes it does but it points to the same object so we should probably just get rid of the hpx attribute.

So I would suggest to support passing in a SparseArray, and to also add SparseArray in all so that it shows up in the docs.

I think the current implementation of __init__ should already support passing in a SparseArray but the docstring should probably be updated to reflect that. I've also expanded the SparseAarray class documentation a bit so it probably makes sense to have it in __all__ now.

But also a factory method on the maps, or even keeping it in init would be OK for me.

I'd prefer having a single method to create an empty map object. Creating an empty map is a very common operation and requiring two function calls seems unnecessarily complex. In practice MapBase.create is a factory method that already partially serves this purpose -- however it doesn't accept a geometry argument but creates a geometry using the method arguments. Maybe we could move the logic in __init__ to a from_geom method that takes a geometry object and a data type.

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Oct 29, 2017

Member

This should be OK now, I've managed to get the BasicFermiLATImageMaker tests to pass again, results are almost identical to before. I'm merging this PR now so that it doesn't block further work on gammapy.maps. For the image maker - if there is a bug, please report it and I'll try to fix, although I would probably try to rewrite it directly into clean cube computation and then integration over energy if possible, instead of spending more time on that image maker.

There are a few PSF-related tests failing with Astropy dev:
https://travis-ci.org/gammapy/gammapy/jobs/294595354#L2597
But I'm pretty sure those fails are unrelated to this PR (because they aren't showing up in any other build).

Member

cdeil commented Oct 29, 2017

This should be OK now, I've managed to get the BasicFermiLATImageMaker tests to pass again, results are almost identical to before. I'm merging this PR now so that it doesn't block further work on gammapy.maps. For the image maker - if there is a bug, please report it and I'll try to fix, although I would probably try to rewrite it directly into clean cube computation and then integration over energy if possible, instead of spending more time on that image maker.

There are a few PSF-related tests failing with Astropy dev:
https://travis-ci.org/gammapy/gammapy/jobs/294595354#L2597
But I'm pretty sure those fails are unrelated to this PR (because they aren't showing up in any other build).

@cdeil cdeil merged commit c421ba9 into gammapy:master Oct 29, 2017

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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment