Skip to content

Commit

Permalink
Fix npred cube tutorial
Browse files Browse the repository at this point in the history
  • Loading branch information
adonath committed Nov 3, 2016
1 parent 2e59980 commit 08790fa
Show file tree
Hide file tree
Showing 3 changed files with 16 additions and 16 deletions.
10 changes: 3 additions & 7 deletions docs/tutorials/npred/npred_general.py
Expand Up @@ -24,7 +24,7 @@ def prepare_images():
exposure_cube.data = Quantity(exposure_cube.data.value, 'cm2 s')

# Re-project background cube
repro_bg_cube = background_model.reproject_to(exposure_cube)
repro_bg_cube = background_model.reproject(exposure_cube)

# Define energy band required for output
energies = EnergyBounds([10, 500], 'GeV')
Expand All @@ -38,12 +38,8 @@ def prepare_images():
offset_max=Angle(3, 'deg'))

# Counts data
counts_data = fits.open(counts_file)[0].data
counts_wcs = WCS(fits.open(counts_file)[0].header)
counts_cube = SkyCube(data=Quantity(counts_data, ''),
wcs=counts_wcs,
energy=energies)
counts_cube = counts_cube.reproject_to(npred_cube, projection_type='nearest-neighbor')
counts_cube = SkyCube.read(counts_file, format='fermi-counts')
counts_cube = counts_cube.reproject(npred_cube)

counts = counts_cube.data[0]
model = convolved_npred_cube.data[0]
Expand Down
12 changes: 9 additions & 3 deletions gammapy/cube/core.py
Expand Up @@ -358,7 +358,7 @@ def spectral(self):
data = np.nansum(np.nansum(self.data, axis=1), axis=1)
return CountsSpectrum(data=data, energy=self.energy)

def lookup(self, position, energy, interpolation=None):
def lookup(self, position, energy, interpolation=False):
"""Differential flux.
Parameters
Expand All @@ -379,8 +379,14 @@ def lookup(self, position, energy, interpolation=None):

z, y, x = self.wcs_skycoord_to_pixel(position, energy)

return self.data[np.rint(z).astype('int'), np.rint(y).astype('int'),
np.rint(x).astype('int')]
if interpolation:
shape = z.shape
pix_coords = np.column_stack([x.flat, y.flat, z.flat])
vals = self._interpolate(pix_coords)
return vals.reshape(shape)
else:
return self.data[np.rint(z).astype('int'), np.rint(y).astype('int'),
np.rint(x).astype('int')]

def show(self, viewer='mpl', ds9options=None, **kwargs):
"""
Expand Down
10 changes: 4 additions & 6 deletions gammapy/cube/utils.py
Expand Up @@ -45,19 +45,17 @@ def compute_npred_cube(flux_cube, exposure_cube, energy_bins,

energy_centers = energy_bins.log_centers
wcs = exposure_cube.wcs
coordinates = exposure_cube.coordinates()
lon = coordinates.data.lon
lat = coordinates.data.lat
coordinates = exposure_cube.spatial.coordinates()
solid_angle = exposure_cube.spatial.solid_angle()

solid_angle = exposure_cube.solid_angle
npred_cube = np.zeros((len(energy_bins) - 1,
exposure_cube.data.shape[1], exposure_cube.data.shape[2]))
for i in range(len(energy_bins) - 1):
energy_bound = energy_bins[i:i + 2]
int_flux = flux_cube.integral_flux_image(energy_bound, integral_resolution)
int_flux = Quantity(int_flux.data, '1 / (cm2 s sr)')
exposure = Quantity(exposure_cube.flux(lon, lat,
energy_centers[i]).value, 'cm2 s')
energy = energy_centers[i] * np.ones(coordinates.shape)
exposure = Quantity(exposure_cube.lookup(coordinates, energy, interpolation=True), 'cm2 s')
npred_image = int_flux * exposure * solid_angle
npred_cube[i] = npred_image.to('')
npred_cube = np.nan_to_num(npred_cube)
Expand Down

0 comments on commit 08790fa

Please sign in to comment.