Add spectrum and flux points to HGPS catalog

gammapy/catalog/fermi.py

@@ -201,17 +201,11 @@ def fetch_fermi_extended_sources(catalog):
 
 
 class SourceCatalogObject3FGL(SourceCatalogObject):
-    """One source from the Fermi-LAT 3FGL catalog.
     """
-
-    x_bins_edges = Quantity([30, 100, 300, 1000, 3000, 10000, 100000], 'MeV')
-
-    x_bins = Quantity(x_bins_edges, 'MeV')
-
-    x_cens = EnergyBounds(x_bins).log_centers
-
-    y_labels = ['Flux30_100', 'Flux100_300', 'Flux300_1000',
-                'Flux1000_3000', 'Flux3000_10000', 'Flux10000_100000']
+    One source from the Fermi-LAT 3FGL catalog.
+    """
+    _ebounds = EnergyBounds(Quantity([100, 300, 1000, 3000, 10000, 100000], 'MeV'))
+    _ebounds_suffix = ['100_300', '300_1000', '1000_3000', '3000_10000', '10000_100000'] 
 
     def __str__(self):
         """Print default summary info string"""
@@ -236,6 +230,57 @@ def summary(self):
 
         return ss
 
+    @property
+    def spectrum(self):
+        raise NotImplementedError
+
+    @property
+    def flux_points_differential(self):
+        """
+        Get `~gammapy.spectrum.DifferentialFluxPoints` for a 3FGL source
+        """
+        from ..spectrum import DifferentialFluxPoints
+
+        energy = self._ebounds.log_centers
+
+        nuFnu = self._get_flux_values('nuFnu', 'erg cm-2 s-1')
+        diff_flux = (nuFnu * energy ** -2).to('erg-1 cm-2 s-1')
+
+        # Get relativ error on integral fluxes
+        int_flux_points = self.flux_points_integral
+        diff_flux_err_hi = diff_flux * int_flux_points['INT_FLUX_ERR_HI_%'] / 100
+        diff_flux_err_lo = diff_flux * int_flux_points['INT_FLUX_ERR_LO_%'] / 100
+
+        return DifferentialFluxPoints.from_arrays(energy=energy, diff_flux=diff_flux,
+                                                  diff_flux_err_lo=diff_flux_err_lo,
+                                                  diff_flux_err_hi=diff_flux_err_hi)
+
+    @property
+    def flux_points_integral(self):
+        """
+        Get `~gammapy.spectrum.IntegralFluxPoints` for a 3FGL source
+
+        Parameters
+        ----------
+        source : dict
+            3FGL source
+        """
+        from ..spectrum import IntegralFluxPoints
+
+        flux = self._get_flux_values()
+        flux_err = self._get_flux_values('Unc_Flux')
+
+        return IntegralFluxPoints.from_arrays(self._ebounds, flux, flux + flux_err[:, 1],
+                                              flux + flux_err[:, 0])
+
+    def _get_flux_values(self, prefix='Flux', unit='cm-2 s-1'):
+        if prefix not in ['Flux', 'Unc_Flux', 'nuFnu']:
+            raise ValueError("Must be one of the following: 'Flux', 'Unc_Flux', 'nuFnu'")
+
+        values = [self.data[prefix + _] for _ in self._ebounds_suffix]
+        return Quantity(values, unit)
+
+
     def plot_lightcurve(self, ax=None):
         """Plot lightcurve.
         """
@@ -244,6 +289,7 @@ def plot_lightcurve(self, ax=None):
         ax = plot_fermi_3fgl_light_curve(self.name, ax=ax)
         return ax
 
+
     def plot_spectrum(self, ax=None):
         """Plot spectrum.
         """

gammapy/catalog/hess.py

@@ -12,6 +12,7 @@
 from __future__ import absolute_import, division, print_function, unicode_literals
 import os
 from collections import OrderedDict
+from astropy.units import Quantity
 from astropy.io import fits
 from astropy.table import Table
 from astropy.coordinates import Angle
@@ -299,6 +300,49 @@ def _summary_associations(self):
         ss += 'List of associated objects: {}\n'.format(associations)
         return ss
 
+    @property
+    def spectrum(self):
+        """
+        `SpectralFitResult` object.
+        """
+        from ..spectrum import SpectrumFitResult
+
+        model = 'PowerLaw'
+        parameters = {}
+        parameters['index'] = Quantity(self.data['Index_Spec_PL'], '')
+        parameters['reference'] = Quantity(self.data['Energy_Spec_PL_Pivot'], 'TeV')
+        parameters['norm'] = Quantity(self.data['Flux_Spec_PL_Diff_Pivot'],
+                                      's^-1 cm^-2 TeV^-1')
+
+        parameter_errors = {}
+        parameter_errors['index'] = Quantity(self.data['Index_Spec_PL_Err'], '')
+        parameter_errors['norm'] = Quantity(self.data['Flux_Spec_PL_Diff_Pivot_Err'],
+                                            's^-1 cm^-2 TeV^-1')
+        parameter_errors['reference'] = Quantity(0, 'TeV')
+
+        ebounds = Quantity([self.data['Energy_Range_Spec_Lo'],
+                            self.data['Energy_Range_Spec_Hi']], 'TeV')
+
+        return SpectrumFitResult(model, parameters, parameter_errors,
+                                 fit_range=ebounds)
+
+    @property
+    def flux_points(self):
+        """
+        Flux points.
+        """
+        from ..spectrum import DifferentialFluxPoints
+
+        energy = Quantity(self.data['Flux_Points_Energy'], 'TeV')
+        diff_flux = Quantity(self.data['Flux_Points_Flux'], 's^-1 cm^-2 TeV^-1')
+        energy_err_hi = Quantity(self.data['Flux_Points_Energy_Err_Hi'], 'TeV')
+        energy_err_lo = Quantity(self.data['Flux_Points_Energy_Err_Lo'], 'TeV')
+        diff_flux_err_hi = Quantity(self.data['Flux_Points_Flux_Err_Hi'], 's^-1 cm^-2 TeV^-1')
+        diff_flux_err_lo = Quantity(self.data['Flux_Points_Flux_Err_Lo'], 's^-1 cm^-2 TeV^-1')
+        return DifferentialFluxPoints.from_arrays(energy, diff_flux, energy_err_hi,
+                                                  energy_err_lo, diff_flux_err_hi,
+                                                  diff_flux_err_lo)
+
 
 class SourceCatalogHGPS(SourceCatalog):
     """HESS Galactic plane survey (HGPS) source catalog.

gammapy/catalog/tests/test_hess.py

@@ -58,3 +58,12 @@ def test_str(self):
         ss = self.source.__str__()
         assert 'Source name          : HESS J1825-137' in ss
         assert 'Component HGPSC 065:' in ss
+
+    def test_spectrum(self):
+        source = self.cat['HESS J1825-137']
+        reference = ('Fit result info \n--------------- \nModel: PowerLaw'
+                     ' \nParameters: \n\t index     : 2.38 +/- 0.03 \n\t norm'
+                     '      : (17.17 +/- 0.57) x 1e-12 / (cm2 s TeV)\n\t reference'
+                     ' : 1.16 +/- 0.00 TeV\n')
+
+        assert str(source.spectrum) == reference

gammapy/spectrum/flux_point.py

@@ -76,38 +76,6 @@ def from_arrays(cls, energy, diff_flux, energy_err_hi=None,
         t['DIFF_FLUX_ERR_LO'] = diff_flux_err_lo
         return cls(t)
 
-    @classmethod
-    def from_3fgl(cls, source, x_method='log_center'):
-        """Get `~gammapy.spectrum.DifferentialFluxPoints` for a 3FGL source
-
-        Parameters
-        ----------
-        sourcename : dict
-            3FGL source
-        """
-        ebounds = EnergyBounds([100, 300, 1000, 3000, 10000, 100000], 'MeV')
-
-        if x_method == 'log_center':
-            energy = ebounds.log_centers
-        else:
-            raise ValueError('Undefined x method: {}'.format(x_method))
-        fluxkeys = ['nuFnu100_300', 'nuFnu300_1000', 'nuFnu1000_3000', 'nuFnu3000_10000', 'nuFnu10000_100000']
-        temp_fluxes = [source.data[_] for _ in fluxkeys]
-        energy_fluxes = Quantity(temp_fluxes, 'erg cm-2 s-1')
-        diff_fluxes = (energy_fluxes * energy ** -2).to('erg-1 cm-2 s-1')
-
-        # Get relativ error on integral fluxes
-        int_flux_points = IntegralFluxPoints.from_3fgl(source)
-        val = int_flux_points['INT_FLUX'].quantity
-        rel_error_hi = int_flux_points['INT_FLUX_ERR_HI'] / val
-        rel_error_lo = int_flux_points['INT_FLUX_ERR_LO'] / val
-
-        diff_fluxes_err_hi = diff_fluxes * rel_error_hi
-        diff_fluxes_err_lo = diff_fluxes * rel_error_lo
-
-        return cls.from_arrays(energy=energy, diff_flux=diff_fluxes,
-                               diff_flux_err_lo=diff_fluxes_err_lo,
-                               diff_flux_err_hi=diff_fluxes_err_hi)
 
     def plot(self, ax=None, energy_unit='TeV',
              flux_unit='cm-2 s-1 TeV-1', energy_power=0, **kwargs):
@@ -140,10 +108,16 @@ def plot(self, ax=None, energy_unit='TeV',
         yh = self['DIFF_FLUX_ERR_HI'].quantity.to(flux_unit).value
         yl = self['DIFF_FLUX_ERR_LO'].quantity.to(flux_unit).value
         y, yh, yl = np.asarray([y, yh, yl]) * np.power(x, energy_power)
+
+        xh = self['ENERGY_ERR_HI'].quantity.to(energy_unit).value
+        xl = self['ENERGY_ERR_LO'].quantity.to(energy_unit).value
+
         flux_unit = Unit(flux_unit) * np.power(Unit(energy_unit), energy_power)
-        ax.errorbar(x, y, yerr=(yl, yh), **kwargs)
+        ax.errorbar(x, y, yerr=(yl, yh), xerr=(xl, xh), **kwargs)
         ax.set_xlabel('Energy [{}]'.format(energy_unit))
         ax.set_ylabel('Flux [{}]'.format(flux_unit))
+        ax.set_xscale("log", nonposx='clip')
+        ax.set_yscale("log", nonposy='clip')
         return ax
 
 
@@ -174,32 +148,10 @@ def from_arrays(cls, ebounds, int_flux, int_flux_err_hi=None,
         t['INT_FLUX'] = int_flux
         t['INT_FLUX_ERR_HI'] = int_flux_err_hi
         t['INT_FLUX_ERR_LO'] = int_flux_err_lo
-        return cls(t)
-
-    @classmethod
-    def from_3fgl(cls, source):
-        """Get `~gammapy.spectrum.IntegralFluxPoints` for a 3FGL source
 
-        Parameters
-        ----------
-        source : dict
-            3FGL source
-        """
-        ebounds = EnergyBounds([100, 300, 1000, 3000, 10000, 100000], 'MeV')
-        fluxkeys = ['Flux100_300', 'Flux300_1000', 'Flux1000_3000', 'Flux3000_10000', 'Flux10000_100000']
-        temp_fluxes = [source.data[_] for _ in fluxkeys]
-
-        fluxerrkeys = ['Unc_Flux100_300', 'Unc_Flux300_1000', 'Unc_Flux1000_3000', 'Unc_Flux3000_10000', 'Unc_Flux10000_100000']
-
-        temp_fluxes_err_hi = [source.data[_][1] for _ in fluxerrkeys]
-        temp_fluxes_err_lo = [-1 * source.data[_][0] for _ in fluxerrkeys]
-
-        int_fluxes = Quantity(temp_fluxes, 'cm-2 s-1')
-        int_fluxes_err_hi = Quantity(temp_fluxes_err_hi, 'cm-2 s-1')
-        int_fluxes_err_lo = Quantity(temp_fluxes_err_lo, 'cm-2 s-1')
-
-        return cls.from_arrays(ebounds, int_fluxes, int_fluxes_err_hi,
-                               int_fluxes_err_lo)
+        t['INT_FLUX_ERR_HI_%'] = 100 * int_flux_err_hi / int_flux
+        t['INT_FLUX_ERR_LO_%'] = 100 * int_flux_err_lo / int_flux
+        return cls(t)
 
     @classmethod
     def from_2fhl(cls, source):

gammapy/spectrum/results.py

@@ -506,6 +506,8 @@ def plot(self, ax=None, energy_unit='TeV', energy_range=None,
         ax.plot(x, y, **kwargs)
         ax.set_xlabel('Energy [{}]'.format(energy_unit))
         ax.set_ylabel('Flux [{}]'.format(flux_unit))
+        ax.set_xscale("log", nonposx='clip')
+        ax.set_yscale("log", nonposy='clip')
         return ax
 
     def plot_butterfly(self, ax=None, energy_unit='TeV', energy_range=None,
@@ -571,6 +573,33 @@ def _get_x(self, energy_range, n_points):
 
         return x
 
+    def __str__(self):
+        """
+        Summary info string.
+        """
+        info = 'Fit result info \n'
+        info += '--------------- \n'
+        info += 'Model: {} \n'.format(self.spectral_model)
+        info += 'Parameters: \n'
+
+        for name in sorted(self.parameters):
+            val = self.parameters[name]
+            err = self.parameter_errors[name]
+            _ = dict(name=name, val=val, err=err)
+            if name == 'norm':
+                _['val'] = _['val'].to('1e-12 cm^-2 TeV^-1 s^-1')
+                _['err'] = _['err'].to('1e-12 TeV^-1 cm^-2 s^-1')
+                info += '\t {name:10s}: ({val.value:.2f} +/- {err.value:.2f}) x {val.unit}\n'.format(**_)
+            else:
+                info += '\t {name:10s}: {val.value:.2f} +/- {err.value:.2f} {val.unit}\n'.format(**_)
+        return info
+
+    def info(self):
+        """
+        Print summary info.
+        """
+        print(str(self))
+
 
 class SpectrumStats(Result):
     """Class summarizing basic spectral parameters

gammapy/spectrum/tests/test_flux_point.py

@@ -194,7 +194,7 @@ def test_3fgl_flux_points():
     source = cat_3fgl['3FGL J0018.9-8152']
     index = source.data['Spectral_Index']
 
-    diff_points = DifferentialFluxPoints.from_3fgl(source)
+    diff_points = source.flux_points_differential
     fluxes = diff_points['DIFF_FLUX'].quantity
     energies = diff_points['ENERGY'].quantity
     eflux = fluxes * energies ** 2
@@ -205,7 +205,7 @@ def test_3fgl_flux_points():
     assert_allclose(actual, diff_desired, rtol=1e-4)
 
 
-    int_points = IntegralFluxPoints.from_3fgl(source)
+    int_points = source.flux_points_integral
     actual = int_points['INT_FLUX'].quantity.to('cm-2 s-1').value
 
     int_desired = [9.45680e-09, 1.94538e-09, 4.00020e-10, 1.26891e-10, 7.24820e-11]