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 SkyImageList #582

Merged
merged 4 commits into from Jul 1, 2016

Conversation

Projects
None yet
3 participants
@OlgaVorokh
Member

OlgaVorokh commented Jun 19, 2016

Initial pull request for SkyImageList -- bridge between SkyMap and SkyCube
Now it's possible to do

cube = SkyCube.read()
sky_image_list = cube.to_image_list()
cube2 = sky_image_list.to_cube()

 cube2 == cube1

and to access individual skymaps like 'sky_image_list.skymaps[i]`

Hi @cdeil,
please take a look and say how close it to what you have imagined and what other functionality should I add :)

def test_to_cube():
sky_cube_original = FermiGalacticCenter.diffuse_model()
image_list = sky_cube_original.to_image_list()
assert np.array_equal(image_list.skymaps[0].data, sky_cube_original.data[0])

This comment has been minimized.

@cdeil

cdeil Jun 19, 2016

Member

Is it possible to just use this here?

assert image_list.skymaps[0].data == sky_cube_original.data[0]

If no, could you please call numpy.testing.assert_array_equal instead of assert np.array_equal? (that should give a better error message in case it's not the same).

This comment has been minimized.

@OlgaVorokh

OlgaVorokh Jun 19, 2016

Member

Comparing arrays with == gives
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all(), thanks for pointing me assert_array_equal!

This comment has been minimized.

@cdeil

cdeil Jun 19, 2016

Member

Also: please move this assert down a little.

Tests are easier to read if they follow the pattern:

# setup
# execution
# asserts

instead of mixing execution and asserts.

@cdeil cdeil added the feature label Jun 19, 2016

@cdeil cdeil added this to the 0.5 milestone Jun 19, 2016

@cdeil cdeil self-assigned this Jun 19, 2016

@cdeil

This comment has been minimized.

Member

cdeil commented Jun 19, 2016

@OlgaVorokh - Thanks!

I will have to think about this and discuss with @adonath before giving comments.
The main question is what the use cases are for this class, compared to SkyMapCollection (related: #578) and whether those should be merged or not.

For now, the main use case for this class will be to implement PSF convolution of a cube.

from gammapy.datasets import FermiGalacticCenter
import numpy as np

This comment has been minimized.

@cdeil

cdeil Jun 19, 2016

Member

You will have to add requires_data here so that the test is skipped on builds that don't have access to test data.

assert np.array_equal(image_list.skymaps[0].data, sky_cube_original.data[0])
sky_cube_restored = image_list.to_cube()
for key in set(sky_cube_original.__dict__.keys() + sky_cube_restored.__dict__.keys()):

This comment has been minimized.

@cdeil

cdeil Jun 19, 2016

Member

I think this code that's looping over keys in the data __dict__ objects is too complicated.
By reading the test it's not clear to me what the keys will be and what exactly happens.

I would prefer to have simple assert statements on the individual data members of the cube
(see https://gammapy.readthedocs.io/en/latest/api/gammapy.cube.SkyCube.html), checking that it's the same.
I.e. roughtly:

# assert cube1.data == cube2.data
# assert cube1.wcs == cube2.wcs
# assert cube1.energy == cube2.energy

(I'm not sure where you can really just check with == or where a more complex assert statement is needed.)

def to_image_list(self):
""" Writes SkyCube to `gammapy.image.SkyImageList` as a list of `gammapy.image.SkyMap` objects.
"""
skymaps = [SkyMap(self.name, Quantity(slice_data, self.data.unit), self.wcs) for slice_data in self.data]

This comment has been minimized.

@cdeil

cdeil Jun 19, 2016

Member

What's the wcs of the resulting skymaps like here?
It should be two-dimensional.
But I think it's three-dimensional because you just copy over the WCS from the cube?

In any case, can you please split this out into a method

cube.sky_image(slice_pos)

that returns an image?
(and add a test for that method that checks that the image WCS is 2-dimensional)

This comment has been minimized.

@OlgaVorokh

OlgaVorokh Jun 19, 2016

Member

Good point!

I printed the wcs that I got from FermiGalacticCenter.diffuse_model() (immediately, for skycube, before converting it) and got

WCS Keywords
Number of WCS axes: 2
CTYPE : 'GLON-CAR'  'GLAT-CAR'  
CRVAL : 0.0  0.0  
CRPIX : 31.5  11.5  
PC1_1 PC1_2  : 1.0  0.0  
PC2_1 PC2_2  : 0.0  1.0  
CDELT : 0.5  0.5  
NAXIS    : 61 21

Is this an error in FermiGalacticCenter.diffuse_model() data files, should I look into this? The code is


    @staticmethod
    def diffuse_model():
        """Diffuse model (`~gammapy.data.SkyCube`)"""
        from ..cube import SkyCube
        filename = FermiGalacticCenter.filenames()['diffuse_model']
        return SkyCube.read(filename)

So as far as I understood from tutorial, the SkyCube usually uses celestial wcs -- longtitute, latitude and cubeface. Can I rely that cubeface is the first axis or should I check that explicitly in cube.sky_image(slice_pos)?

This comment has been minimized.

@cdeil

cdeil Jun 20, 2016

Member

Actually I now see that the SkyCube.wcs is already the 2-dimensional celestial WCS slice ... this is done in SkyCube.read:
https://github.com/gammapy/gammapy/blob/master/gammapy/cube/core.py#L158

I don't know if I coded that up a year or two ago or if someone changed it to be like this in the meantime.

In any case the class docstring is correct ... it should say that the wcs member is the 2-dim celestial WCS here:
https://github.com/gammapy/gammapy/blob/master/gammapy/cube/core.py#L45

@OlgaVorokh So you can just pass the cube WCS to the image WCS, and don't have to do any WCS computation or axis order check. Could you please update the cube docstring to a correct one-sentence statement what WCS is in this PR?

@@ -324,6 +320,22 @@ def solid_angle_image(self):
sky_map = SkyMap.from_image_hdu(image_hdu)
return sky_map.solid_angle()
def sky_image(self, slicepos):
"""Slice 3-dim cube into a 2-dim `~gammapy.image.SkyMap along slicepos`.

This comment has been minimized.

@cdeil

cdeil Jun 21, 2016

Member

The closing backtick is in the wrong place here ... it should be:

Slice out a 2-dim `~gammapy.image.SkyMap` at a given energy slice position.

One important question is whether this method should create a "view" or a "copy".

A view is more efficient, but this can also be very surprising to users if they create an image this way, modify it, and then the underlying cube changes.
I think the default behaviour should be the "safe version" of copy=True to create a decoupled image, and then there should be an option copy=False that one can use for the few cases where the copy is too costly and for experts.
@adonath - Do you agree?

Also -- the current behaviour of this method is not clear to me. I'm pretty sure name and wcs are by reference (i.e. coupled, not a copy). For data, I think self.data[slicepos] is a view and
This will produce a "view", i.e. create a SkyMap where the name, data, unit and Quantity(self.data[slicepos], self.data.unit) also doesn't make a copy ... but here I'm not sure.

@OlgaVorokh - Could you play with this (or an equivalent little test case with a Quantity in IPython) a little until you understand when a copy is made and when it isn't, and then add the copy=False option to this method and implement it?
I don't think the behaviour needs to be verified via tests ... this is something we can do via code review. But if you want to add a test that establishes that copy=False and copy=True work as expected, that's fine too.

This comment has been minimized.

@cdeil

cdeil Jun 21, 2016

Member

Actually, @OlgaVorokh - Could you please split out the SkyCube.sky_image method into a separate PR and do this before continuing with this one?

In the same PR you can also remove the standalone cube_to_image function

def cube_to_image(cube, slicepos=None):

and update the one place where it's called: SkyCube.solid_angle.

This comment has been minimized.

@adonath

adonath Jun 21, 2016

Member

For SkyMap, there is a .copy() method, maybe it would make sense to have that for SkyCube as well? Then one could copy the cube before creating the SkyMap. But there should definitely be the option copy=False, that gives you a reference ('view') to the data. Which is useful, if you want to access large data with the SkyImageList interface instead of SkyCube.

This comment has been minimized.

@OlgaVorokh

OlgaVorokh Jun 22, 2016

Member

@cdeil Would it be okay to call SkyMap.copy() on a newly created copy instead of working with each parameter and determining if it's mutable?
I think the overhead is negligible as SkyMap.init consist of just several "=" operators

@@ -324,6 +320,10 @@ def solid_angle_image(self):
sky_map = SkyMap.from_image_hdu(image_hdu)
return sky_map.solid_angle()
def sky_image(self, slicepos):

This comment has been minimized.

@adonath

adonath Jun 21, 2016

Member

Wouldn't it be more natural to give an energy (e.g. center energy) here instead of of slicepos?

This comment has been minimized.

@cdeil

cdeil Jun 21, 2016

Member

The version with the index is definitely needed.

E.g. when computing the solid angle image, one would simply do sky_image(slicepos=0).solid_angle() and when computing the PSF-convolved cube one would iterate over energy slice indices.

But yes, the use case of getting the slice for a given energy is also valid.
It could be added here by putting energy_idx=None, energy=None and then checking that the caller supplied exactly one.
Or we could rename SkyCube.flux to SkyCube.evaluate and have a "neareast neighbor" option that doesn't do interpolation.

Maybe we table that question for now and get the basic version in first?

@OlgaVorokh OlgaVorokh referenced this pull request Jun 24, 2016

Merged

Add SkyCube.sky_image() #586

@OlgaVorokh OlgaVorokh force-pushed the OlgaVorokh:skyimagelist branch from 4563afd to 8862a53 Jun 28, 2016

@@ -324,6 +320,24 @@ def solid_angle_image(self):
sky_map = SkyMap.from_image_hdu(image_hdu)
return sky_map.solid_angle()
def sky_image(self, idx_energy, copy=True):
"""Slice 3-dim cube into a 2-dim `~gammapy.image.SkyMap` along slicepos.

This comment has been minimized.

@cdeil

cdeil Jun 29, 2016

Member

The SkyCube.sky_image method was merged in #586, I can now see it in master:
https://github.com/gammapy/gammapy/blob/master/gammapy/cube/core.py#L324

I don't know why it still shows up in the diff for this PR on Github!?

@OlgaVorokh Can you please rebase this PR and this will refresh the diff?

@OlgaVorokh OlgaVorokh force-pushed the OlgaVorokh:skyimagelist branch from 8862a53 to 30e1b6c Jun 29, 2016

@OlgaVorokh

This comment has been minimized.

Member

OlgaVorokh commented Jun 29, 2016

@cdeil, I have rebased this branch on master

from . import catalog
from ..image import SkyMap
from ..irf import EnergyDependentTablePSF
from ..cube import SkyCube

This comment has been minimized.

@cdeil

cdeil Jun 30, 2016

Member

@OlgaVorokh - Can you please put this SkyCube import as a delayed import inside the SkyImageList.to_cube method?
Recently we discussed a Gammapy sub-package dependency structure, i.e. which sub-packages are lower-level and which are higher-level, to avoid circular imports. (I still have to write that up.)
gammapy.image is lower than gammapy.cube, i.e. should mostly avoid importing things from gammapy.cube, and where needed, do it as a delayed import.

@@ -0,0 +1,44 @@
from __future__ import absolute_import, division, print_function, unicode_literals

This comment has been minimized.

@cdeil

cdeil Jun 30, 2016

Member

Please add the license comment line that's at the top of any file in the Gammapy to gammapy/image/lists.py and gammapy/image/tests/test_lists.py.

from astropy.units import Quantity
from astropy.wcs import WCS
from ..utils.testing import requires_dependency, requires_data
from . import catalog

This comment has been minimized.

@cdeil

cdeil Jun 30, 2016

Member

Some of these imports are unused, no?
E.g. these can be removed:

from . import catalog
from ..irf import EnergyDependentTablePSF

If you edit in PyCharm this will be pointed out to you directly.
Alternatively, you can run e.g. https://pypi.python.org/pypi/pyflakes on your files from the command line, and it will tell you which imports are unused.

@cdeil

This comment has been minimized.

Member

cdeil commented Jun 30, 2016

I've left some minor comments inline.
Once those are addressed, I would suggest we merge this.

It's still not clear how this SkyImageList and the existing SkyMapCollection are related. Probably they should be merged into one class. But figuring this out will take longer and is already under discussion in #578.

@adonath - Are you OK if we merge this for now and then continue the image list discussion / implementation in other issues / pull requests?

@adonath

This comment has been minimized.

Member

adonath commented Jun 30, 2016

@cdeil For me the SkyImageList class will be the main class to handle images of one data type e.g. counts, exposure etc. with different binning for every energy band (or even the same energy band, as a kind of "scale-space" image list). This is not reflected in the API yet (e.g. the to_cube() method probably needs an additional call to .reproject() to reproject all the images to a common wcs).
Whereas the SkyMapCollection is meant to be a generic higher level data container for maps, cubes or image lists of different type. I.e. bundle counts, background, exposure data for simple FITS I/O and convenient bunch style attribute access for users. And to bundle the data in a single container for input and output of other functions such as compute_ts_map or adaptive_smooth.

Yes, I'm totally fine with merging this PR and further develop the API of SkyImageList in future PRs.

@OlgaVorokh

This comment has been minimized.

Member

OlgaVorokh commented Jun 30, 2016

@cdeil I fixed all you asked me

@cdeil

This comment has been minimized.

Member

cdeil commented Jul 1, 2016

Looks good to me. Thanks!

@cdeil cdeil merged commit 02842b1 into gammapy:master Jul 1, 2016

2 checks passed

continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment