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 code to make model images from a source catalog #160

Merged
merged 6 commits into from Aug 22, 2014

Conversation

@ellisowen
Copy link
Contributor

@ellisowen ellisowen commented Aug 4, 2014

Incomplete, work in progress - I made a pull request to allow discussion about the structure & implementation.

Note that some functions also use stuff used in npred_cube2 pull request (which is why I rebased against it in this PR).

To do (Ellis):

  • Check output type & format looks sensible
  • Ensure output values are sensible (Flux normalization)
  • Write unit tests
  • Add doc strings
  • Add short tutorial
@cdeil
Copy link
Member

@cdeil cdeil commented Aug 4, 2014

Let's try and finish #151 first ... reviewing a PR that contains commits from another PR is difficult, because the diff is large and it's not clear what should be reviewed in this PR.

One thing you can do now is make format the task list above as described here so that you can tick the items off as you complete them.

from astropy.coordinates import Angle


def reference_image(resolution, latitudes, longitudes, center=[0,0], units='ph/cm2/s/sr'):

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014
Member

You should not add a new reference_image function here.
Use (and improve if necessary) the existing make_header and / or make_empty_image functions.

energy = Quantity([10, 500], 'GeV')):
"""TODO
"""
reference = reference_image(resolution, lat_range, lon_range, center, units='ph/cm2/s/sr') #Check these units!

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014
Member

You should pass in the reference image instead of passing all the parameters needed to construct it.

I'm never sure what's best, to pass an ImageHDU or a Header or a WCS and shape ... see the reproject function as an example that can handle all of those.
You don't need to support that here ... just pick one option ... maybe ImageHDU?

reference = reference_image(resolution, lat_range, lon_range, center, units='ph/cm2/s/sr') #Check these units!
lons, lats = coordinates(reference)
wcs = WCS(reference.header)
total_point_image = GammaSpectralCube(data = Quantity(np.array([reference.data]), ''), wcs = wcs, energy = energy)

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014
Member

The usual stuff ... add a few empty lines ... break lines or use temp variables to have line length <~ 80 characters to make the code more readable.

table.meta['Energy Bins'] = energy
return table

if __name__ == '__main__':

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014
Member

Remove.

return out_cube


def catalog_table(catalog, ebands='No'):

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014
Member

So you're just adding things to the Fermi catalog table?
Then maybe make this a private function and optionally call it at the end of fetch_fermi_catalog?

@cdeil
Copy link
Member

@cdeil cdeil commented Aug 4, 2014

Yes, adding this to the existing npred tutorial or making a new one that creates npred images for the Fermi catalog sources would be nice (and show that this actually works)

@ellisowen
Copy link
Contributor Author

@ellisowen ellisowen commented Aug 7, 2014

Positioning and generating images works OK, but flux normalization is off by about 30 orders of magnitude between the point sources and diffuse sources... I still need to understand why.

@ellisowen
Copy link
Contributor Author

@ellisowen ellisowen commented Aug 7, 2014

To save time, would it be OK if I just implement the point source stuff for one energy band? This is all that's needed for the proceeding, and I think figuring out what the problem is with the flux of the extended sources, and extending the code to work with cubes would be a lot of work...

@cdeil
Copy link
Member

@cdeil cdeil commented Aug 7, 2014

Yes, one energy band is OK.

I've been thinking lately that probably using lists of 2D images as much as possible could be better than using 3D cubes. Then one could simply concatenate lists to e.g. work on Fermi and HESS data simultaneously.

@ellisowen
Copy link
Contributor Author

@ellisowen ellisowen commented Aug 9, 2014

This now works + tests + tutorial added.

The Travis CI build should pass once npred_cube2 and vela_corrections pull requests have been merged (at which point I will also be able to clear up the commit history here).

There are a few semi-implemented options and functions included which might be useful to implement in the future but would take too long to do now. I left these with 'NotImplementedErrors' for now and TODO comments.

@cdeil
Copy link
Member

@cdeil cdeil commented Aug 20, 2014

Now that #151 is merged into master, could you please rebase this PR?

@ellisowen ellisowen force-pushed the ellisowen:image_from_catalog branch 2 times, most recently from 08afaa4 to f1230ba Aug 20, 2014
@ellisowen
Copy link
Contributor Author

@ellisowen ellisowen commented Aug 20, 2014

OK, rebased & ready for review... although I see there is new functionality catalog to cube in gammapy.astro.population - I'm not sure if this duplicates what I implement here?

image = catalog_image(reference, psf, catalog='simulation', source_type='point',
total_flux=True, sim_table=table)

# Plot

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

Please change the plotting for sky images to use alpy so that the images have coordinate axes.

Here's an example that works in this case:

# Plot
from aplpy import FITSFigure
fig = FITSFigure(image.to_fits()[0])
fig.show_grayscale(stretch='sqrt', interpolation='none')
fig.add_colorbar()
@cdeil
Copy link
Member

@cdeil cdeil commented Aug 21, 2014

You have to add

from .catalog import *

to gammapy/image/__init__.py.
Then please check that the API docs look OK by running python setup.py build_sphinx.

from gammapy.astro.population import add_observed_parameters
from gammapy.datasets import FermiGalacticCenter
from gammapy.image import make_empty_image
from gammapy.image.catalog import catalog_image

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

You should use:

from gammapy.image import catalog_image
@ellisowen ellisowen force-pushed the ellisowen:image_from_catalog branch from f1230ba to aa500f6 Aug 21, 2014
@ellisowen
Copy link
Contributor Author

@ellisowen ellisowen commented Aug 21, 2014

I wasn't able to check the aplpy images because of the issues I mentioned before, but I guess if they didn't work Travis CI would fail?

Otherwise, this should now address your comments.

@cdeil
Copy link
Member

@cdeil cdeil commented Aug 21, 2014

I've built the html docs locally and put the online for discussion / reference.
The tutorial is here and the two functions added here are catalog_image and catalog_table.

For the tutorial, could you maybe make a 300 x 100 pix image so that one sees more of the Galactic plane?
I don't remember what aplpy issue you had, but you need to resolve that anyways to make nice plots for your report and the proceeding, no? Email me if you can't get it to work.

I'll leave some inline comments ... one more round of code review.

from aplpy import FITSFigure
from gammapy.datasets import FermiGalacticCenter
from gammapy.image import make_empty_image
from gammapy.image.catalog import catalog_image

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

Import from gammapy.image, not from gammapy.image.catalog

from gammapy.image.catalog import catalog_image
from gammapy.irf import EnergyDependentTablePSF

# Define image size

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

It's good if example code is short, so I'd suggest to delete these variables write here and simply write:

reference = make_empty_image(nxpix=100, nypix=100, binsz=1)
Catalog & Simulation Images
===========================

The `~gammapy.image.catalog.catalog_image` method allows the production of

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

`~gammapy.image.catalog.catalog_image`

should be changed to

`~gammapy.image.catalog_image`
from astropy.wcs import WCS
from astropy.units import Quantity
from astropy.table import Table
from ..datasets.load import fetch_fermi_extended_sources, fetch_fermi_catalog

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

Import from ..dataset, not from ..datasets.load here.
(always use the location that's shown in the the API html docs)

source_spec_cube = GammaSpectralCube(data=Quantity(np.array([source.data]), ''),
wcs=source_wcs, energy=energy)
new_source_cube = source_spec_cube.reproject_to(reference_cube)
# TODO: Fix this hack

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

Can you resolve this TODO?
(or explain what it's about and I'll try to help)

This comment has been minimized.

@ellisowen

ellisowen Aug 21, 2014
Author Contributor

The flux values are about 30 orders of magnitude larger than I expected for the extended image, and I don't know why...
However, I don't need this function anymore so I don't think there's any point in implementing it now - shall I put a not implemented error or just remove the code?

This comment has been minimized.

@cdeil

cdeil Aug 22, 2014
Member

At some point we do want to add the extended sources ... so let's leave this as-is and figure it out later.

x, y = wcs.wcs_world2pix(lon, lat, 0)
xi, yi = x.astype(int), y.astype(int)
new_image[yi, xi] = new_image[yi, xi] + flux
if total_flux == True:

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

Remove == True ... it's not needed.

new_image[yi, xi] = new_image[yi, xi] + flux
if total_flux == True:
factor = source_table['flux'].sum() / new_image.sum()
return new_image * factor, energies

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

Move the return statement at the end.
Instead modify or create new variables in the if statement and always return at the end.

Functions that have multiple return paths are always more complex and if it can easily be avoided you should.

Energy dependent Table PSF object for image convolution.
catalog : str
Flag which source catalog is to be used to create the image.
Options: {'1FHL', '2FGL', 'simulation'}

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

The docstring convention for string options is

catalog : {'1FHL', '2FGL', 'simulation'}

and not to put

catalog : str

and the list the options in the description.

There's a few places where this needs to be changed ... please read over the diff on Github and fix this everywhere.

point_source = _source_image(catalog, reference_cube, total_flux=True)[0]
new_image = extended + point_source
else:
raise ValueError

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

Add empty newline after such an if block ... better readable.



def catalog_table(catalog, ebands=False):
""" Creates catalog table from published source catalog.

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

Remove empty space """ Creates ... -> """Creates ...

return out_cube


def catalog_table(catalog, ebands=False):

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

ebands -> in_energy_bands


def catalog_table(catalog, ebands=False):
""" Creates catalog table from published source catalog.

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

From the function name and docstring alone it's not possible to figure out what this function really does.
Add a sentence about what it does or when it should be used.

glon = cat_table['GLON'][source]
glat = cat_table['GLAT'][source]
# Different from here between each catalog because of different catalog header names
if catalog == ('1FHL' or 'simulation'):

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

I don't think this makes sense.

In [6]: catalog = 'simulation'

In [7]: catalog == ('1FHL' or 'simulation')
Out[7]: False

Maybe you want this?

if catalog in ['1FHL', 'simulation']:

Is this code covered by a unit test?

This comment has been minimized.

@ellisowen

ellisowen Aug 21, 2014
Author Contributor

Yes - although I don't explicitly test the 'simulation' option as this does exactly the same as '1FHL'. I only include both options to make it a bit clearer for the user.

# Different from here between each catalog because of different catalog header names
if catalog == ('1FHL' or 'simulation'):
energy = Quantity([10, 30, 100, 500], 'GeV')
if ebands == False:

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

Use this instead:

if not ebands:

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

Or better yet, use

if ebands:

and switch the if and else clause.

Flux_1000_3000 = cat_table['Flux1000_3000'][source]
Flux_3000_10000 = cat_table['Flux3000_10000'][source]
Flux_10000_100000 = cat_table['Flux10000_100000'][source]
row = dict(Source_Type='PointSource', Source_Name=source_name,

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

I find such things more readable if there's one entry per line (even if it takes more lines).

Flux_100_300=Flux_100_300, Flux_300_1000=Flux_300_1000,
Flux_1000_3000=Flux_1000_3000, Flux_3000_10000=Flux_3000_10000,
Flux_10000_100000=Flux_10000_100000)
data.append(row)

This comment has been minimized.

@cdeil

cdeil Aug 21, 2014
Member

Add empty line here for better readability.

@ellisowen ellisowen force-pushed the ellisowen:image_from_catalog branch from aa500f6 to ead8a93 Aug 21, 2014
@ellisowen
Copy link
Contributor Author

@ellisowen ellisowen commented Aug 21, 2014

OK, I think I got all your comments.
I left one or two questions in line.

The Travis CI failures don't currently make sense to me; the errors are quite cryptic (& don't occur when I run setup.py locally)

(Oh, and I've got aplpy to work on my laptop, so I'll probably make a PR in the next day or so to change the npred cube tutorial images to use it)

@ellisowen ellisowen force-pushed the ellisowen:image_from_catalog branch from ead8a93 to 5ce3eac Aug 21, 2014
@cdeil
Copy link
Member

@cdeil cdeil commented Aug 22, 2014

I also don't know what this travis-ci failure is about:
https://travis-ci.org/gammapy/gammapy/jobs/33231332
I've re-started all travis-ci builds for this repo ... lets' see if the issue disappears.

@cdeil
Copy link
Member

@cdeil cdeil commented Aug 22, 2014

Travis-ci passed ... merging this now.

Thanks!

cdeil added a commit that referenced this pull request Aug 22, 2014
Image from catalog
@cdeil cdeil merged commit c4b48ed into gammapy:master Aug 22, 2014
1 check passed
1 check passed
continuous-integration/travis-ci The Travis CI build passed
Details
@cdeil cdeil changed the title Image from catalog Add code to make model images from a source catalog 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
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

None yet

2 participants
You can’t perform that action at this time.