Merge pull request #3480 from adonath/flux_profile_estimator

Implement FluxProfileEstimator

gammapy/datasets/core.py

@@ -287,6 +287,29 @@ def slice_by_energy(self, energy_min, energy_max):
 
         return self.__class__(datasets=datasets)
 
+    def to_spectrum_datasets(self, region):
+        """Extract spectrum datasets for the given region.
+
+        Parameters
+        ----------
+        region : `~regions.SkyRegion`
+            Region definition.
+
+        Returns
+        --------
+        datasets : `Datasets`
+            List of `~gammapy.datasets.SpectrumDataset`
+        """
+        datasets = Datasets()
+
+        for dataset in self:
+            spectrum_dataset = dataset.to_spectrum_dataset(
+                on_region=region, name=dataset.name
+            )
+            datasets.append(spectrum_dataset)
+
+        return datasets
+
     @property
     # TODO: make this a method to support different methods?
     def energy_ranges(self):

gammapy/datasets/map.py

@@ -1329,9 +1329,9 @@ def to_spectrum_dataset(self, on_region, containment_correction=False, name=None
                 )
                 dataset.exposure.quantity *= containment.reshape(geom.data_shape)
 
-        kwargs = {}
+        kwargs = {"name": name}
 
-        for name in [
+        for key in [
             "counts",
             "edisp",
             "mask_safe",
@@ -1340,12 +1340,12 @@ def to_spectrum_dataset(self, on_region, containment_correction=False, name=None
             "gti",
             "meta_table",
         ]:
-            kwargs[name] = getattr(dataset, name)
+            kwargs[key] = getattr(dataset, key)
 
         if self.stat_type == "cash":
             kwargs["background"] = dataset.background
 
-        return SpectrumDataset(**kwargs)
+        return  SpectrumDataset(**kwargs)
 
     def to_region_map_dataset(self, region, name=None):
         """Integrate the map dataset in a given region.
@@ -2350,9 +2350,12 @@ def to_spectrum_dataset(self, on_region, containment_correction=False, name=None
         """
         from .spectrum import SpectrumDatasetOnOff
 
-        dataset = super().to_spectrum_dataset(on_region, containment_correction, name)
+        dataset = super().to_spectrum_dataset(
+            on_region=on_region, containment_correction=containment_correction, name=name
+        )
+
+        kwargs = {"name": name}
 
-        kwargs = {}
         if self.counts_off is not None:
             kwargs["counts_off"] = self.counts_off.get_spectrum(
                 on_region, np.sum, weights=self.mask_safe

gammapy/estimators/excess_map.py

@@ -151,24 +151,12 @@ def run(self, dataset):
         Parameters
         ----------
         dataset : `~gammapy.datasets.MapDataset` or `~gammapy.datasets.MapDatasetOnOff`
-            input dataset
+            Map dataset
 
         Returns
         -------
-        images : dict
-            Dictionary containing result correlated maps. Keys are:
-
-                * counts : correlated counts map
-                * background : correlated background map
-                * excess : correlated excess map
-                * ts : TS map
-                * sqrt_ts : sqrt(delta TS), or Li-Ma significance map
-                * err : symmetric error map (from covariance)
-                * flux : flux map. An exposure map must be present in the dataset to compute flux map
-                * errn : negative error map
-                * errp : positive error map
-                * ul : upper limit map
-
+        maps : `FluxMaps`
+            Flux maps
         """
         if not isinstance(dataset, MapDataset):
             raise ValueError("Unsupported dataset type. Excess map is not applicable to 1D datasets.")

gammapy/estimators/flux.py

@@ -167,5 +167,8 @@ def run(self, datasets):
         models[self.source].spectral_model = model
         datasets.models = models
         result.update(super().run(datasets, model.norm))
+
+        # TODO: find a cleaner way of inlcuding the npred_excess info
+        datasets.models[self.source].spectral_model.norm.value = result["norm"]
         result.update(self.estimate_npred_excess(datasets=datasets))
         return result

gammapy/estimators/flux_map.py

@@ -700,13 +700,19 @@ def from_stack(cls, maps, axis, meta=None):
         """
         reference = maps[0]
         data = {}
+
         for quantity in reference.available_quantities:
             data[quantity] = Map.from_stack([_._data[quantity] for _ in maps], axis=axis)
 
         if meta is None:
             meta = reference.meta.copy()
 
-        gti = GTI.from_stack([_.gti for _ in maps])
+        gtis = [_.gti for _ in maps if _.gti]
+
+        if gtis:
+            gti = GTI.from_stack(gtis)
+        else:
+            gti = None
 
         return cls(
             data=data,

gammapy/estimators/flux_point.py

@@ -385,6 +385,26 @@ def to_table(self, sed_type="likelihood", format="gadf-sed", formatted=False):
                 data = getattr(self, quantity, None)
                 if data:
                     table[quantity] = data.quantity.squeeze()
+        elif format == "profile":
+            x_axis = self.geom.axes["projected-distance"]
+
+            tables = []
+            for idx, (x_min, x_max) in enumerate(x_axis.iter_by_edges):
+                table_flat = Table()
+                table_flat["x_min"] = [x_min]
+                table_flat["x_max"] = [x_max]
+                table_flat["x_ref"] = [(x_max + x_min) / 2]
+
+                fp = self.slice_by_idx(slices={"projected-distance": idx})
+                table = fp.to_table(sed_type=sed_type, format="gadf-sed")
+
+                for column in table.columns:
+                    table_flat[column] = table[column][np.newaxis]
+
+                tables.append(table_flat)
+
+            table = vstack(tables)
+
         else:
             raise ValueError(f"Not a supported format {format}")
 
@@ -637,17 +657,15 @@ def run(self, datasets):
 
         Parameters
         ----------
-        datasets : list of `~gammapy.datasets.Dataset`
+        datasets : `~gammapy.datasets.Datasets`
             Datasets
 
         Returns
         -------
         flux_points : `FluxPoints`
             Estimated flux points.
         """
-        # TODO: remove copy here...
-        datasets = Datasets(datasets).copy()
-
+        datasets = Datasets(datasets=datasets)
         rows = []
 
         for energy_min, energy_max in progress_bar(
@@ -692,6 +710,5 @@ def estimate_flux_point(self, datasets, energy_min, energy_max):
         datasets_sliced = datasets.slice_by_energy(
             energy_min=energy_min, energy_max=energy_max
         )
-
         datasets_sliced.models = datasets.models.copy()
         return super().run(datasets=datasets_sliced)

gammapy/estimators/profile.py

@@ -6,9 +6,139 @@
 from astropy.convolution import Box1DKernel, Gaussian1DKernel
 from astropy.coordinates import Angle
 from astropy.table import Table
+from regions import CircleAnnulusSkyRegion
+from gammapy.datasets import Datasets
+from gammapy.maps import MapAxis
+from gammapy.modeling.models import PowerLawSpectralModel, SkyModel
+from .flux_point import FluxPoints, FluxPointsEstimator
 from .core import Estimator
 
-__all__ = ["ImageProfile", "ImageProfileEstimator"]
+
+__all__ = ["ImageProfile", "ImageProfileEstimator", "FluxProfileEstimator"]
+
+
+class FluxProfileEstimator(FluxPointsEstimator):
+    """Estimate flux profiles
+
+    Parameters
+    ----------
+    regions : list of `SkyRegion`
+        regions to use
+    spectrum : `~gammapy.modeling.models.SpectralModel` (optional)
+        Spectral model to compute the fluxes or brightness.
+        Default is power-law with spectral index of 2.
+    **kwargs : dict
+        Keywords forwarded to the `FluxPointsEstimator` (see documentation
+        there for further description of valid keywords)
+
+    Examples
+    --------
+    This example shows how to compute a counts profile for the Fermi galactic
+    center region::
+
+        from astropy import units as u
+        from astropy.coordinates import SkyCoord
+        from gammapy.data import GTI
+        from gammapy.estimators import FluxProfileEstimator
+        from gammapy.utils.regions import make_orthogonal_rectangle_sky_regions
+        from gammapy.datasets import MapDataset
+        from gammapy.maps import RegionGeom
+
+        # load example data
+        dataset = MapDataset.read("$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc.fits.gz", name="fermi_dataset")
+
+        # configuration
+        dataset.gti = GTI.create("0s", "1e7s", "2010-01-01")
+
+        # creation of the boxes and axis
+        start_pos = SkyCoord("-1d", "0d", frame='galactic')
+        end_pos = SkyCoord("1d", "0d", frame='galactic')
+
+        regions = make_orthogonal_rectangle_sky_regions(
+            start_pos=start_pos,
+            end_pos=end_pos,
+            wcs=dataset.counts.geom.wcs,
+            height=2.*u.deg,
+            nbin=21
+        )
+
+        # set up profile estimator and run
+        prof_maker = FluxProfileEstimator(
+            regions=regions, energy_edges=[10, 2000] * u.GeV,
+        )
+        fermi_prof = prof_maker.run(dataset)
+
+        print(fermi_prof)
+
+    """
+    tag = "FluxProfileEstimator"
+
+    def __init__(self, regions, spectrum=None, **kwargs):
+        if len(regions) <= 1:
+            raise ValueError("Please provide at least two regions for flux profile estimation.")
+
+        self.regions = regions
+
+        if spectrum is None:
+            spectrum = PowerLawSpectralModel()
+
+        self.spectrum = spectrum
+        super().__init__(**kwargs)
+
+    @property
+    def projected_distance_axis(self):
+        """Get projected distance from the first region.
+
+        For normal region this is defined as the distance from the
+        center of the region. For annulus shaped regions it is the
+        mean between the inner and outer radius.
+
+        Returns
+        -------
+        axis : `MapAxis`
+            Projected distance axis
+        """
+        distances = []
+        center = self.regions[0].center
+
+        for idx, region in enumerate(self.regions):
+            if isinstance(region, CircleAnnulusSkyRegion):
+                distance = (region.inner_radius + region.outer_radius) / 2.0
+            else:
+                distance = center.separation(region.center)
+
+            distances.append(distance)
+
+        return MapAxis.from_nodes(
+            u.Quantity(distances, "deg"), name="projected-distance"
+        )
+
+    def run(self, datasets):
+        """Run flux profile estimation
+
+        Parameters
+        ----------
+        datasets : list of `~gammapy.datasets.MapDataset`
+            Map datasets.
+
+        Returns
+        -------
+        profile : `~gammapy.estimators.FluxPoints`
+            Profile flux points.
+        """
+        datasets = Datasets(datasets=datasets)
+
+        maps = []
+
+        for region in self.regions:
+            datasets_to_fit = datasets.to_spectrum_datasets(region=region)
+            datasets_to_fit.models = SkyModel(self.spectrum, name="test-source")
+            fp = super().run(datasets_to_fit)
+            maps.append(fp)
+
+        return FluxPoints.from_stack(
+            maps=maps, axis=self.projected_distance_axis,
+        )
 
 
 # TODO: implement measuring profile along arbitrary directions