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 SED from Cube function #166

Merged
merged 5 commits into from Aug 22, 2014

Conversation

Projects
None yet
2 participants
@ellisowen
Contributor

ellisowen commented Aug 15, 2014

Made a pull request while I work on this in case you have comments.

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Aug 17, 2014

Ok, this is now ready - please review.
The Travis CI failures don't appear to be related to anything I've changed here.

@cdeil cdeil referenced this pull request Aug 20, 2014

Merged

Add image profile function #167

else:
errors = np.zeros_like(values)
if flux_type == 'Differential':

This comment has been minimized.

@cdeil

cdeil Aug 20, 2014

Member

Always use lower-case for string options: 'Differential' -> 'differential'

return table
else:
raise ValueError

This comment has been minimized.

@cdeil

cdeil Aug 20, 2014

Member

Give the usual message to the ValueError (search the code for other examples if it's not clear what I mean).

Quantity(values, cube.data.unit),
Quantity(errors * values, cube.data.unit)],
names=('ENERGY', 'DIFF_FLUX', 'DIFF_FLUX_ERR'))
return table

This comment has been minimized.

@cdeil

cdeil Aug 20, 2014

Member

Put the return table statement from here and from the next elif clause at the end of the function outside the if statement.

@ellisowen ellisowen force-pushed the ellisowen:sed branch from b25d2b9 to 35788dd Aug 20, 2014

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Aug 20, 2014

OK, this should address your comments

as floats. A 1x2 array specifying the maximum and minimum latitudes to
include in determining the SED points, which should be within the
spatial bounds of the cube.
flux_type : String, {'differential', 'integral'}

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014

Member

Remove "String, " here ... just the part in the braces is the right format.

include in determining the SED points, which should be within the
spatial bounds of the cube.
flux_type : String, {'differential', 'integral'}
Specify whether input cube includes Differential or Integral fluxes.

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014

Member

Differential -> differential
Integral -> integral

counts : `~gammapy.spectral_cube.GammaSpectralCube`
Counts cube to allow Poisson errors to be calculated. If not provided,
a standard_error should be provided, or zero errors will be returned.
(Optional).

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014

Member

The "Optional"` should be removed from the description and be put in this line at the end after a comma:

counts :  `~gammapy.spectral_cube.GammaSpectralCube`, optional
values.append(bin.value)
values = np.array(values)
if errors == True:

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014

Member

Remove the == True part.

values = np.array(values)
if errors == True:
if counts == None:

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014

Member

Use this instead:

if counts is None:

(please turn on PEP8 checking in your Eclipse editor ... it will show you this as a warning, so you'll get it right automatically)

int_flux_err=errors * values)
else:
raise ValueError('Unknown selection: {0}'.format(flux_type))

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014

Member

selection -> flux_type

@cdeil

This comment has been minimized.

Member

cdeil commented Aug 21, 2014

Can you please add a tutorial where you make a GLON, GLAT profile and SED for the Fermi Galactic center example?
(maybe in a new PR if you prefer)

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Aug 21, 2014

  • Can you please add a tutorial where you make a GLON, GLAT profile and SED for the Fermi Galactic center example?
    (maybe in a new PR if you prefer)

Sure, I will add this in a new PR (when I am producing the plots for the report anyway so can just copy my code across when it's ready).

@ellisowen ellisowen force-pushed the ellisowen:sed branch from 35788dd to dd0940a Aug 21, 2014

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Aug 21, 2014

OK, this should now address your comments.

@@ -7,7 +7,7 @@
from .diffuse import *
from .fitter import *
from . import fitting_utils
from .flux_point import *
from .flux_point import compute_differential_flux_points

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014

Member

Why did you change this?
compute_differential_flux_points is listed in __all__ in flux_point.py.

We do not want to have a duplication of the manual __all__ lists in the __init__.py files ... that would be hard to maintain (and there's no reason to do it).

This comment has been minimized.

@ellisowen

ellisowen Aug 21, 2014

Contributor

Ah, this was previously to fix some circular import errors. I resolved the problem but forgot to switch this back. Now done.

Returns
-------
table : `~astropy.table.Table`
A Spectral Energy table of Energies, Differential Fluxes and

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014

Member

No need to capitalize the first letter ... use this:

A spectral energy table of energies, differential fluxes and errors.
def cube_sed(cube, lats, lons, flux_type='differential', counts=None, mask=None,
errors=False, standard_error=0.1, spectral_index=2.3):
"""Creates SED from GammaSpectralCube within given GLAT and GLON range.

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014

Member

GLAT -> lat, GLON -> lon.

Passing only a mask that defines an arbitrary region would be more flexible, i.e. remote the lats and lons parameters.
If it doesn't exist yet, you can add a utility function to gammapy.image to make a mask for a rectangle in lon and lat.

if flux_type == 'differential':
energy = cube.energy
table = Table([energy,

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014

Member

I find it more readable if you do

table = Table()
table['ENERGY'] = energy
table['DIFF_FLUX'] = Quantity(values, cube.data.unit)
table['DIFF_FLUX_ERR'] = Quantity(errors * values, cube.data.unit)

It's one simple statement per line instead of a complex statement across four lines,
and it has column name and value together.

counts = FermiGalacticCenter.diffuse_model()
counts.data = np.ones_like(counts.data)
mask = np.zeros_like(spec_cube.data[0]).value

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014

Member

Please add a comment here ... what kind of mask does this create?
Does this include all pixels?
Having some pixels included and some not would make a better test, no?

mask = np.zeros_like(spec_cube.data[0]).value
lats, lons = spec_cube.spatial_coordinate_images
lat = [lats.min(), lats.max()]

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014

Member

This is brittle code ... I'm not sure if this takes the whole image or if the edge pixels are removed.

If you change the cube_sed function following my suggestion above of only taking a mask as input,
you could simply remove this code if you want the full image.

And if you want a rectangular sub-part, you call a utility function that makes this mask.
I don't know if it exists, but it could be

lon_lat_rectangle_mask(lons, lats, lon_min=None, lon_max=None, lat_min=None, lat_max=None)

where the default None means no limit.

mask = np.zeros_like(spec_cube.data[0]).value
lats, lons = spec_cube.spatial_coordinate_images

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014

Member

I think the order is incorrect here ... according to the docs GammaSpectralCube.spatial_coordinate_images returns lon, lat and not lat, lon as you write here.
See https://gammapy.readthedocs.org/en/latest/api/gammapy.spectral_cube.GammaSpectralCube.html?#gammapy.spectral_cube.GammaSpectralCube.spatial_coordinate_images

Please check and fix GammaSpectralCube.spatial_coordinate_images if you have reason to think it might be incorrect.

This comment has been minimized.

@ellisowen

ellisowen Aug 22, 2014

Contributor

Yes, this was my mistake - GammaSpectralCube.spatial_coordinate_images is OK.

sed_table3 = cube_sed(spec_cube, lat, lon, flux_type='differential',
errors=True, standard_error = 0.1)
assert_allclose(sed_table3['DIFF_FLUX'].data, 12810 * np.ones(30))

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014

Member

The * np.ones(30) here should be surperfluous, because numpy broadcasting takes care of this, no?

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014

Member

There's a few other lines nearby where np.ones can simply be removed, I think.

This comment has been minimized.

@ellisowen

ellisowen Aug 22, 2014

Contributor

Yes, now changed this!

@cdeil

This comment has been minimized.

Member

cdeil commented Aug 21, 2014

I've left another bunch of inline comments, the most important ones are that lons and lats should be dropped as arguments to this function (all you need is a mask) and that you switched lon and lat in the unit test, possibly making the reference results incorrect (unless you were selecting all pixels in the image all the time, which wouldn't be a good test).

@ellisowen ellisowen force-pushed the ellisowen:sed branch 2 times, most recently from 97fc695 to 91a554b Aug 22, 2014

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Aug 22, 2014

OK, how is this?

lats : array_like
Array of galactic latitude values.
lon_min : float, optional
Minimum galactic longitude of rectangular mask.

This comment has been minimized.

@cdeil

cdeil Aug 22, 2014

Member

Remove 'galactic' ... this can be used for any coordinate system.

Parameters
----------
lons : array_like
Array of galactic longitude values.

This comment has been minimized.

@cdeil

cdeil Aug 22, 2014

Member

Again ... remove "galactic" everywhere.

@@ -195,3 +196,11 @@ def test_wcs_histogram2d():
assert measure.lookup(image, 0, 0, world=False) == 1 + 3
assert measure.lookup(image, 1, 0, world=False) == 2
def test_lon_lat_rectangle_mask():

This comment has been minimized.

@cdeil

cdeil Aug 22, 2014

Member

Add tests for the default None cases for the limits also.
You can use this command to make sure every line in your function is tested:

python setup.py test -V --coverage -t gammapy/image/tests/test_utils.py
Returns
-------
mask : array_like

This comment has been minimized.

@cdeil

cdeil Aug 22, 2014

Member

array_like ->

`numpy.ndarray`
Parameters
----------
lons : array_like

This comment has been minimized.

@cdeil

cdeil Aug 22, 2014

Member

I don't think this will work for array_like inputs like e.g. a list of lists, because you call numpy methods on it.
So change to numpy.array here for the type of the inputs.

@cdeil

This comment has been minimized.

Member

cdeil commented Aug 22, 2014

This is mostly OK ... I left some more inline comments.

@ellisowen ellisowen force-pushed the ellisowen:sed branch from 91a554b to b28557e Aug 22, 2014

@ellisowen ellisowen force-pushed the ellisowen:sed branch from b28557e to a8c3c11 Aug 22, 2014

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Aug 22, 2014

I think this addresses your comments

@cdeil

This comment has been minimized.

Member

cdeil commented Aug 22, 2014

Thanks!

cdeil added a commit that referenced this pull request Aug 22, 2014

Merge pull request #166 from ellisowen/sed
SED from Cube code

@cdeil cdeil merged commit 5022770 into gammapy:master Aug 22, 2014

1 check passed

continuous-integration/travis-ci The Travis CI build passed
Details

@cdeil cdeil changed the title from SED from Cube code to Add SED from Cube function Apr 8, 2015

@cdeil cdeil added the feature label Apr 8, 2015

@cdeil cdeil added this to the 0.1 milestone Apr 8, 2015

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment