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 SkyEllipse model #2046

Merged

Conversation

Projects
2 participants
@luca-giunti
Copy link
Contributor

luca-giunti commented Feb 21, 2019

Added elliptical model defined and normalized on the sphere. The doc includes an example that produces an illustrative plot
ellipse

@adonath adonath self-requested a review Feb 22, 2019

@adonath adonath self-assigned this Feb 22, 2019

@adonath adonath added the feature label Feb 22, 2019

@adonath adonath added this to In progress in Modeling via automation Feb 22, 2019

@adonath adonath added this to the 0.11 milestone Feb 22, 2019

@adonath
Copy link
Member

adonath left a comment

Thanks @luca-giunti! I've left a few in-line comments mainly concerning the code style. I really like the example plot. It's very illustrative. We should consider adding one for the other models as well.

My main comment would be to add a second test, using a larger ellipse radius to better check the normalization.

:math:`\text{lon}_0`: `lon` coordinate for the center of the ellipse.
lat_0 : `~astropy.coordinates.Latitude`
:math:`\text{lat}_0`: `lat` coordinate for the center of the ellipse.
semi_major : `~astropy.coordinates.Angle`

This comment has been minimized.

Copy link
@adonath

adonath Feb 22, 2019

Member

Maybe it's an option to rename it to r_0 to be consistent with SkyDisk model, for e = 0...I have no strong opinion here.

This comment has been minimized.

Copy link
@luca-giunti

luca-giunti Feb 22, 2019

Author Contributor

Thank you for all your comments! So yes, using r_0 might be an option. I think it might lead to some confusion, however. Indeed, if the model is wrongly passed the minor semiaxis it misbehaves

Eccentricity of the ellipse (:math:`0< e< 1`).
theta : `~astropy.coordinates.Angle`
:math:`\theta`:
Rotation angle of the major semiaxis. The

This comment has been minimized.

Copy link
@adonath

adonath Feb 22, 2019

Member

Here is a bad line break...

Examples
--------

This comment has been minimized.

Copy link
@adonath

adonath Feb 22, 2019

Member

The title and underline have to be aligned to be correctly rendered by sphinx...

:include-source:
import numpy as np
from gammapy.image.models.core import *

This comment has been minimized.

Copy link
@adonath

adonath Feb 22, 2019

Member

Please replace by the explicit import from gammapy.image.models import SkyEllipse. In Python it's actually never recommended to use the * import, because it can silently overwrite variable and function names. Imagine you would use from math import * and from numpy import *, this would lead to a completely messy namespace...

This comment has been minimized.

Copy link
@luca-giunti

luca-giunti Feb 22, 2019

Author Contributor

Thank you! very instructive

from gammapy.maps import Map, WcsGeom
import matplotlib.pyplot as plt
from matplotlib import rc
rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})

This comment has been minimized.

Copy link
@adonath

adonath Feb 22, 2019

Member

This modifies the global rc state of matplotlib, that means if a user does a second plot in the same Python session, those changes will also apply to the second one, which might not be desired. Instead you could use the rc_context context manager. You'll find examples on the linked page, as well as in the Gammapy code e.g. https://github.com/gammapy/gammapy/blob/master/gammapy/maps/base.py#L998

]
)

def find_foci(lon_0, lat_0, semi_major, e, theta):

This comment has been minimized.

Copy link
@adonath

adonath Feb 22, 2019

Member

I think it's not super-useful to have this as a separate method, just because I don't think it will be called by users. It could be in the .evaluate() method as well...do as you prefer.

assert_allclose(
np.sum(mymap.quantity * u.sr ** -1 * m_allsky.solid_angle()), 1, rtol=5.0e-4
)

This comment has been minimized.

Copy link
@adonath

adonath Feb 22, 2019

Member

Maybe add a second test case with an even larger radius for the ellipse (~10 deg) to check that the normalisation is correct?

This comment has been minimized.

Copy link
@luca-giunti

luca-giunti Feb 22, 2019

Author Contributor

Yes sure! I'll do that. Altough I think it may somewhat increase the computation time for the test

This comment has been minimized.

Copy link
@adonath

adonath Feb 22, 2019

Member

Of course, you can choose a coarser pixel size to limit the computation time...

norm = u.Quantity(
SkyEllipse.compute_norm(semi_major, e), unit="sr-1", copy=False
)
result = np.select([in_ellipse], [norm])

This comment has been minimized.

Copy link
@adonath

adonath Feb 22, 2019

Member

I think np.select() does not work properly with Quantity, that's probably why the instance check is needed later. You could simplify the code as following:

norm = SkyEllipse.compute_norm()
return u.Quantity(norm * in_ellipse, "sr-1", copy=False)

This implicitly does the type-cast of the bool array in_ellipse, but it works correctly with Quantity as well.

vals=ell(lon,lat)
mymap=Map.from_geom(m_geom,data=vals.value)
fig,ax,_=mymap.plot()

This comment has been minimized.

Copy link
@adonath

adonath Feb 22, 2019

Member

Please reformat the example according to PEP8 standards (https://www.python.org/dev/peps/pep-0008/). There is probably a plugin avalaible for your editor...

@adonath adonath merged commit dd8b9f9 into gammapy:master Feb 28, 2019

0 of 3 checks passed

Codacy/PR Quality Review Hang in there, Codacy is reviewing your Pull request.
Details
Scrutinizer Running
Details
continuous-integration/travis-ci/pr The Travis CI build is in progress
Details

Modeling automation moved this from In progress to Done Feb 28, 2019

@adonath adonath changed the title Add elliptical model defined on the sphere Add SkyEllipse model Mar 26, 2019

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.