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

Change map smooth widths to match Astropy #1807

Merged
merged 5 commits into from
Sep 20, 2018
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/image/survey_example.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
filename = "$GAMMAPY_EXTRA/datasets/fermi_survey/all.fits.gz"
survey_map = Map.read(filename, hdu="counts")
survey_map.data = survey_map.data.astype("float")
smoothed_map = survey_map.smooth(radius=Angle(0.2, unit="deg"))
smoothed_map = survey_map.smooth(width=Angle(0.1, unit="deg"))

fig = plt.figure(figsize=(15, 8))
xlim = Angle([70, 262], unit="deg")
Expand Down
2 changes: 1 addition & 1 deletion gammapy/data/event_list.py
Original file line number Diff line number Diff line change
Expand Up @@ -659,7 +659,7 @@ def _counts_image(self):
}
m = WcsNDMap.create(**opts)
m.fill_by_coord(self.radec)
m = m.smooth(radius=1)
m = m.smooth(width=1)
return m

def plot_image(self):
Expand Down
4 changes: 2 additions & 2 deletions gammapy/maps/tests/test_wcsnd.py
Original file line number Diff line number Diff line change
Expand Up @@ -489,8 +489,8 @@ def test_convolve_vs_smooth():
m = WcsNDMap.create(binsz=binsz, width=1.05 * u.deg, axes=axes)
m.data[:, :, 10, 10] = 1.

desired = m.smooth(kernel="gauss", radius=0.5 * u.deg, mode="constant")
gauss = Gaussian2DKernel(5).array
desired = m.smooth(kernel="gauss", width=0.5 * u.deg, mode="constant")
gauss = Gaussian2DKernel(10).array
actual = m.convolve(kernel=gauss)
assert_allclose(actual.data, desired.data, rtol=1e-3)

Expand Down
16 changes: 7 additions & 9 deletions gammapy/maps/wcsnd.py
Original file line number Diff line number Diff line change
Expand Up @@ -491,19 +491,20 @@ def _plot_format_allsky(self, ax):
lat.grid(alpha=0.2, linestyle="solid", color="w")
return ax

def smooth(self, radius, kernel="gauss", **kwargs):
def smooth(self, width, kernel="gauss", **kwargs):
"""
Smooth the image (works on a 2D image and returns a copy).

The definition of the smoothing parameter radius is equivalent to the
one that is used in ds9 (see `ds9 smoothing <http://ds9.si.edu/doc/ref/how.html#Smoothing>`_).

Parameters
----------
radius : `~astropy.units.Quantity` or float
width : `~astropy.units.Quantity` or float
Smoothing width given as quantity or float. If a float is given it
interpreted as smoothing width in pixels. If an (angular) quantity
is given it converted to pixels using ``geom.wcs.wcs.cdelt``.
It corresponds to the standard deviation in case of a Gaussian kernel,
the radius in case of a disk kernel, and the side length in case
of a box kernel.
kernel : {'gauss', 'disk', 'box'}
Kernel shape
kwargs : dict
Expand All @@ -518,22 +519,19 @@ def smooth(self, radius, kernel="gauss", **kwargs):
"""
from scipy.ndimage import gaussian_filter, uniform_filter, convolve

if isinstance(radius, u.Quantity):
radius = (radius.to("deg") / self.geom.pixel_scales.mean()).value
if isinstance(width, u.Quantity):
width = (width.to("deg") / self.geom.pixel_scales.mean()).value

smoothed_data = np.empty_like(self.data)

for img, idx in self.iter_by_image():
if kernel == "gauss":
width = radius / 2.
data = gaussian_filter(img, width, **kwargs)
elif kernel == "disk":
width = 2 * radius + 1
disk = Tophat2DKernel(width)
disk.normalize("integral")
data = convolve(img, disk.array, **kwargs)
elif kernel == "box":
width = 2 * radius + 1
data = uniform_filter(img, width, **kwargs)
else:
raise ValueError("Invalid kernel: {!r}".format(kernel))
Expand Down