Repository: gammapy/gammapy
Pull request #5538 (2024-11-04 10:27:07)

Add a helper function to create a global dataset from a list of dataset.
Add a helper function to create a global dataset from a list of dataset.
Follow-up of #5533 that introduces  a function `create_empty_map_dataset_from_irfs`. 

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/gammapy/gammapy/pull/5538?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 94.60%. Comparing base [(`4b03188`)](https://app.codecov.io/gh/gammapy/gammapy/commit/4b0318881c2b4a0d026033840792c6a76404c183?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) to head [(`4cb8d65`)](https://app.codecov.io/gh/gammapy/gammapy/commit/4cb8d659b2d6e2f55794d4d9babe04916ff8251a?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).
> Report is 56 commits behind head on main.

<details><summary>Additional details and impacted files</summary>


```diff
@@           Coverage Diff           @@
##             main    #5538   +/-   ##
=======================================
  Coverage   94.60%   94.60%           
=======================================
  Files         232      232           
  Lines       36293    36350   +57     
=======================================
+ Hits        34335    34390   +55     
- Misses       1958     1960    +2     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/gammapy/gammapy/pull/5538?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).

QRemy: > Maybe you could clarify the scope and use case of the function?
The long term goal is to determine automatically the dataset geoms from the irfs instead of recommending specific defaults in the configuration files of the hli.
It is intended to be used after  `create_empty_map_dataset_from_irfs` (which create one dataset with a binning optimised for each irf)  in case we want a single dataset to stack all observations.

> Does it work in all situations?
Could be refined later if it breaks for special cases.

> Why do you create a dataset in galactic Plate-Carrée projection in the end?
modified, now it uses the frame and projection of the first dataset in the list


Patch:
diff --git a/gammapy/datasets/utils.py b/gammapy/datasets/utils.py
index 15d7fb6d5..a5c88ef27 100644
--- a/gammapy/datasets/utils.py
+++ b/gammapy/datasets/utils.py
@@ -2,7 +2,9 @@
 import numpy as np
 from astropy.coordinates import SkyCoord
 from gammapy.data import Observation
-from gammapy.maps import Map
+import astropy.units as u
+from gammapy.maps import Map, MapAxis, WcsGeom
+from gammapy.maps.utils import _check_width, _check_binsz
 from gammapy.modeling.models.utils import cutout_template_models
 from . import Datasets
 
@@ -287,3 +289,131 @@ def create_map_dataset_from_dl4(data, geom=None, energy_axis_true=None, name=Non
         dataset.edisp.exposure_map.data = dataset.psf.exposure_map.data
 
     return dataset
+
+
+def create_global_dataset(
+    datasets,
+    name=None,
+    position=None,
+    binsz=None,
+    width=None,
+    energy_min=None,
+    energy_max=None,
+    energy_true_min=None,
+    energy_true_max=None,
+    nbin_per_decade=None,
+):
+    """Create an empty dataset encompassing the input datasets.
+
+    Parameters
+    ----------
+    datasets : `~gammapy.datasets.Datasets`
+        Datasets
+    name : str, optional
+        Name of the output dataset. Default is None.
+    position : `~astropy.coordinates.SkyCoord`
+        Center position of the output dataset.
+        Default is None, and the average position is taken.
+    binsz : float or tuple or list, optional
+        Map pixel size in degrees.
+        Default is None, the minimum bin size is taken.
+    width : float or tuple or list or string, optional
+        Width of the map in degrees.
+        Default is None, in that case it is derived as
+        the maximal width + twice the maximal separation between datasets and position.
+    energy_min :  `~astropy.units.Quantity`
+        Energy range.
+        Default is None, the minimum energy is taken.
+    energy_max :  `~astropy.units.Quantity`
+        Energy range.
+        Default is None, the maximum energy is taken.
+    energy_true_min, energy_true_max :  `~astropy.units.Quantity`,  `~astropy.units.Quantity`
+        True energy range. If None, minimum and maximum energies of all geometries is used.
+        Default is None.
+    nbin_per_decade : int
+        number of energy bins per decade.
+        Default is None, the maximum is taken.
+
+    Returns
+    -------
+    datasets : `~gammapy.datasets.MapDatset`
+        Empy global dataset.
+    """
+    from gammapy.datasets import MapDataset
+
+    if position is None:
+        frame = datasets[0].counts.geom.frame
+        positions = SkyCoord(
+            [d.counts.geom.center_skydir.transform_to(frame) for d in datasets]
+        )
+        position = SkyCoord(positions.cartesian.mean(), frame="icrs")
+        position = SkyCoord(position.ra, position.dec, frame="icrs").transform_to(
+            frame
+        )  # drop fake distance
+
+    binsz_list = []
+    width_list = []
+    energy_min_list = []
+    energy_max_list = []
+    energy_true_min_list = []
+    energy_true_max_list = []
+    nbin_per_decade_list = []
+    for d in datasets:
+        binsz_list.append(np.abs(d.counts.geom.pixel_scales).min())
+        width_list.append(
+            d.counts.geom.width.max()
+            + 2 * d.counts.geom.center_skydir.separation(position)
+        )
+        energy_min_list.append(d.counts.geom.axes["energy"].edges.min())
+        energy_max_list.append(d.counts.geom.axes["energy"].edges.max())
+        energy_true_min_list.append(d.exposure.geom.axes["energy_true"].edges.min())
+        energy_true_max_list.append(d.exposure.geom.axes["energy_true"].edges.max())
+        ndecade = np.log10(energy_true_max_list[-1].value) - np.log10(
+            energy_true_min_list[-1].value
+        )
+        nbins = len(d.exposure.geom.axes[0].center)
+        nbin_per_decade_list.append(np.ceil(nbins / ndecade))
+
+    if binsz is None:
+        binsz = np.min(u.Quantity(binsz_list))
+    binsz = _check_binsz(binsz)
+    if width is None:
+        width = np.max(u.Quantity(width_list))
+    width = _check_width(width)
+    energy_true_min = (
+        energy_true_min
+        if energy_true_min is not None
+        else u.Quantity(energy_true_min_list).min()
+    )
+    if energy_true_max is None:
+        energy_true_max = np.max(u.Quantity(energy_true_max_list))
+    if energy_min is None:
+        energy_min = np.min(u.Quantity(energy_min_list))
+    if energy_max is None:
+        energy_max = np.max(u.Quantity(energy_max_list))
+
+    nbin_per_decade = (
+        nbin_per_decade if nbin_per_decade is not None else np.max(nbin_per_decade_list)
+    )
+    energy_axis = MapAxis.from_energy_bounds(
+        energy_min, energy_max, nbin_per_decade, unit="TeV", per_decade=True
+    )
+
+    geom = WcsGeom.create(
+        skydir=position,
+        binsz=binsz,
+        width=width,
+        frame=position.frame,
+        proj=datasets[0].counts.geom.projection,
+        axes=[energy_axis],
+    )
+
+    energy_axis_true = MapAxis.from_energy_bounds(
+        energy_true_min,
+        energy_true_max,
+        nbin_per_decade,
+        unit="TeV",
+        name="energy_true",
+        per_decade=True,
+    )
+    return MapDataset.create(geom=geom, energy_axis_true=energy_axis_true, name=name)


Test patch:
diff --git a/gammapy/datasets/tests/test_utils.py b/gammapy/datasets/tests/test_utils.py
index c81539f63..00932caf5 100644
--- a/gammapy/datasets/tests/test_utils.py
+++ b/gammapy/datasets/tests/test_utils.py
@@ -4,7 +4,7 @@ import numpy as np
 from numpy.testing import assert_allclose
 import astropy.units as u
 from gammapy.datasets import MapDataset
-from gammapy.datasets.utils import apply_edisp, split_dataset
+from gammapy.datasets.utils import apply_edisp, split_dataset, create_global_dataset
 from gammapy.irf import EDispKernel
 from gammapy.maps import Map, MapAxis
 from gammapy.modeling.models import (
@@ -103,3 +103,61 @@ def test_dataset_split():
     assert (~datasets[7].mask_fit.data[0, :, :]).sum() == np.prod(
         pixel_width + 2 * margin_width
     ) - np.prod(pixel_width)
+
+
+@requires_data()
+def test_create_global_dataset():
+    base_dataset = MapDataset.read("$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc.fits.gz")
+    width = 4 * u.deg
+    margin = 1 * u.deg
+    datasets = split_dataset(base_dataset, width, margin)
+
+    global_dataset = create_global_dataset(datasets)
+    assert global_dataset.counts.geom.width.max() > base_dataset.counts.geom.width.max()
+    assert (
+        global_dataset.counts.geom.axes[0].nbin > base_dataset.counts.geom.axes[0].nbin
+    )
+    assert_allclose(
+        global_dataset.counts.geom.axes[0].edges.min(),
+        base_dataset.counts.geom.axes[0].edges.min(),
+    )
+    assert_allclose(
+        global_dataset.counts.geom.axes[0].edges.max(),
+        base_dataset.counts.geom.axes[0].edges.max(),
+    )
+    assert_allclose(
+        global_dataset.exposure.geom.axes[0].edges.min(),
+        base_dataset.exposure.geom.axes[0].edges.min(),
+    )
+    assert_allclose(
+        global_dataset.exposure.geom.axes[0].edges.max(),
+        base_dataset.exposure.geom.axes[0].edges.max(),
+    )
+
+    global_dataset = create_global_dataset(
+        datasets,
+        width=(20, 10) * u.deg,
+        binsz=0.02 * u.deg,
+        name="test",
+        position=datasets[0].counts.geom.center_skydir,
+        energy_min=0.05 * u.TeV,
+        energy_max=0.2 * u.TeV,
+        energy_true_min=0.02 * u.TeV,
+        energy_true_max=0.4 * u.TeV,
+        nbin_per_decade=10,
+    )
+    assert global_dataset.name == "test"
+    assert_allclose(global_dataset.counts.geom.width.max(), 20 * u.deg)
+    assert_allclose(global_dataset.counts.geom.pixel_scales.max(), 0.02 * u.deg)
+    assert_allclose(global_dataset.counts.geom.axes[0].nbin, 7)
+    assert_allclose(
+        global_dataset.counts.geom.center_skydir.separation(
+            datasets[0].counts.geom.center_skydir
+        ).value,
+        0,
+        atol=1e-5,
+    )
+    assert_allclose(global_dataset.counts.geom.axes[0].edges.min(), 0.05 * u.TeV)
+    assert_allclose(global_dataset.counts.geom.axes[0].edges.max(), 0.2 * u.TeV)
+    assert_allclose(global_dataset.exposure.geom.axes[0].edges.min(), 0.02 * u.TeV)
+    assert_allclose(global_dataset.exposure.geom.axes[0].edges.max(), 0.4 * u.TeV)