Implement SkyCube.spectrum() method
@@ -504,8 +504,9 @@ class SkyCube(object):
data = np.zeros_like(self.data[0])
return SkyImage(name=self.name, data=data, wcs=wcs)
- def lookup(self, position, energy, interpolation=False):
- """Differential flux.
+ def lookup(self, position, energy, interpolation=None):
+ """
+ Lookup value in the cube at given sky position and energy.
Parameters
----------
@@ -513,11 +514,13 @@ class SkyCube(object):
Position on the sky.
energy : `~astropy.units.Quantity`
Energy
+ interpolation : {None, 'linear', 'nearest'}
+ Interpolate data values between energies.
Returns
-------
- flux : `~astropy.units.Quantity`
- Differential flux (1 / (cm2 MeV s sr))
+ value : `~astropy.units.Quantity`
+ Value at the given sky position and energy.
"""
# TODO: add interpolation option using NDDataArray
@@ -721,35 +724,60 @@ class SkyCube(object):
images = [self.sky_image(energy) for energy in energies]
return SkyCubeImages(self.name, images, self.wcs, energies)
- def to_spectrum(self, region, weights=None):
+ def spectrum(self, region):
"""
- Integrate spatial dimensions of the cube to give a spectrum.
-
- TODO: give formulas.
+ Extract spectrum in a given sky region.
Parameters
----------
region : `~regions.SkyRegion`
- Region to sum the data in.
- weighting : {None, 'solid_angle'}, optional
- Weighting factor to use.
+ Sky region to extract the spectrum from.
Returns
-------
- spectrum : numpy.array
+ spectrum : `~astropy.table.Table`
Summed spectrum of pixels in the mask.
"""
- raise NotImplementedError
- # TODO: clean up API and implementation and add test
- value = cube.dat
+ spectrum = Table()
- mask = region.contains(sky_image_ref.coordinates())
- if weights == 'solid_angle':
- weights = self.sky_image_ref.solid_angle()
+ # store energy binning
+ energies = self.energies('edges')
+ e_ref = self.energies('center')
+ spectrum['e_min'] = energies[:-1]
+ spectrum['e_max'] = energies[1:]
+ spectrum['e_ref'] = e_ref
- spectrum = np.nansum(data * weights).sum(-1).sum(-1)
+ # mask region and sum
+ mask = self.region_mask(region)
+ value = (self.data * mask.data).sum(-1).sum(-1)
+ spectrum['value'] = value
return spectrum
+ def region_mask(self, region):
+ """
+ Create a boolean cube mask for a region.
+
+ The mask is:
+
+ - ``True`` for pixels inside the region
+ - ``False`` for pixels outside the region
+
+ Parameters
+ ----------
+ region : `~regions.PixelRegion` or `~regions.SkyRegion`
+ Region in pixel or sky coordinates.
+
+ Returns
+ -------
+ mask : `SkyCube`
+ A boolean sky cube mask.
+ """
+ mask = self.sky_image_ref.region_mask(region)
+ data = mask.data * np.ones(self.data.shape, dtype='bool') * u.Unit('')
+ wcs = self.wcs.deepcopy() if self.wcs else None
+ return self.__class__(name=self.name, data=data.astype('bool'), wcs=wcs,
+ energy_axis=self.energy_axis)
+
def write(self, filename, format='fermi-counts', **kwargs):
"""Write to FITS file.
