Merge pull request #4853 from AtreyeeS/dataset_meta

Add MapDatasetMetaData container

gammapy/datasets/__init__.py

@@ -10,6 +10,9 @@
     create_map_dataset_geoms,
 )
 from .simulate import MapDatasetEventSampler, ObservationEventSampler
+from .map import MapDataset, MapDatasetOnOff, create_map_dataset_geoms
+from .metadata import MapDatasetMetaData
+from .simulate import MapDatasetEventSampler
 from .spectrum import SpectrumDataset, SpectrumDatasetOnOff
 
 DATASET_REGISTRY = Registry(
@@ -39,4 +42,5 @@
     "OGIPDatasetReader",
     "SpectrumDataset",
     "SpectrumDatasetOnOff",
+    "MapDatasetMetaData",
 ]

gammapy/datasets/map.py

@@ -25,6 +25,7 @@
 from gammapy.utils.table import hstack_columns
 from .core import Dataset
 from .evaluator import MapEvaluator
+from .metadata import MapDatasetMetaData
 from .utils import get_axes
 
 __all__ = [
@@ -306,6 +307,8 @@ class MapDataset(Dataset):
     meta_table : `~astropy.table.Table`
         Table listing information on observations used to create the dataset.
         One line per observation for stacked datasets.
+    meta : `~gammapy.datasets.MapDatasetMetaData`
+        Associated meta data container
 
 
     Notes
@@ -390,6 +393,7 @@ def __init__(
         gti=None,
         meta_table=None,
         name=None,
+        meta=None,
     ):
         self._name = make_name(name)
         self._evaluators = {}
@@ -420,6 +424,18 @@ def __init__(
         self.gti = gti
         self.models = models
         self.meta_table = meta_table
+        if meta is None:
+            self._meta = MapDatasetMetaData()
+        else:
+            self._meta = meta
+
+    @property
+    def meta(self):
+        return self._meta
+
+    @meta.setter
+    def meta(self, value):
+        self._meta = value
 
     # TODO: keep or remove?
     @property
@@ -1023,6 +1039,9 @@ def stack(self, other, nan_to_num=True):
         elif other.meta_table:
             self.meta_table = other.meta_table.copy()
 
+        if self.meta and other.meta:
+            self.meta.stack(other.meta)
+
     def stat_array(self):
         """Statistic function value per bin given the current model parameters."""
         return cash(n_on=self.counts.data, mu_on=self.npred().data)
@@ -1337,6 +1356,7 @@ def to_hdulist(self):
 
         header = hdu_primary.header
         header["NAME"] = self.name
+        header.update(self.meta.to_header())
 
         hdulist = fits.HDUList([hdu_primary])
         if self.counts is not None:
@@ -1390,6 +1410,7 @@ def from_hdulist(cls, hdulist, name=None, lazy=False, format="gadf"):
         """
         name = make_name(name)
         kwargs = {"name": name}
+        kwargs["meta"] = MapDatasetMetaData.from_header(hdulist["PRIMARY"].header)
 
         if "COUNTS" in hdulist:
             kwargs["counts"] = Map.from_hdulist(hdulist, hdu="counts", format=format)
@@ -2206,6 +2227,8 @@ class MapDatasetOnOff(MapDataset):
         One line per observation for stacked datasets.
     name : str
         Name of the dataset.
+    meta : `~gammapy.datasets.MapDatasetMetaData`
+        Associated meta data container
 
 
     See Also
@@ -2232,6 +2255,7 @@ def __init__(
         mask_safe=None,
         gti=None,
         meta_table=None,
+        meta=None,
     ):
         self._name = make_name(name)
         self._evaluators = {}
@@ -2248,6 +2272,10 @@ def __init__(
         self.models = models
         self.mask_safe = mask_safe
         self.meta_table = meta_table
+        if meta is None:
+            self._meta = MapDatasetMetaData()
+        else:
+            self._meta = meta
 
     def __str__(self):
         str_ = super().__str__()

gammapy/datasets/metadata.py

@@ -0,0 +1,101 @@
+from typing import ClassVar, Literal, Optional, Union
+import numpy as np
+from astropy.coordinates import SkyCoord
+from pydantic import ConfigDict
+from gammapy.utils.metadata import (
+    METADATA_FITS_KEYS,
+    CreatorMetaData,
+    MetaData,
+    ObsInfoMetaData,
+    PointingInfoMetaData,
+)
+
+__all__ = ["MapDatasetMetaData"]
+
+MapDataset_METADATA_FITS_KEYS = {
+    "MapDataset": {
+        "event_types": "EVT_TYPE",
+        "optional": "OPTIONAL",
+    },
+}
+
+METADATA_FITS_KEYS.update(MapDataset_METADATA_FITS_KEYS)
+
+
+class MapDatasetMetaData(MetaData):
+    """Metadata containing information about the Dataset.
+
+    Parameters
+    ----------
+    creation : `~gammapy.utils.CreatorMetaData`, optional
+         The creation metadata.
+    obs_info : list of `~gammapy.utils.ObsInfoMetaData`
+        info about the observation.
+    event_types : list of int or str
+        Event types used in analysis.
+    pointing: list of `~gammapy.utils.PointingInfoMetaData`
+        Telescope pointing directions.
+    optional : dict
+        Additional optional metadata.
+    """
+
+    model_config = ConfigDict(coerce_numbers_to_str=True)
+
+    _tag: ClassVar[Literal["MapDataset"]] = "MapDataset"
+    creation: Optional[CreatorMetaData] = CreatorMetaData()
+    obs_info: Optional[Union[ObsInfoMetaData, list[ObsInfoMetaData]]] = None
+    pointing: Optional[Union[PointingInfoMetaData, list[PointingInfoMetaData]]] = None
+    event_type: Optional[Union[str, list[str]]] = None
+    optional: Optional[dict] = None
+
+    def stack(self, other):
+        kwargs = {}
+        kwargs["creation"] = self.creation
+        return self.__class__(**kwargs)
+
+    @classmethod
+    def _from_meta_table(cls, table):
+        """Create MapDatasetMetaData from MapDataset.meta_table
+
+        Parameters
+        ----------
+            table: `~astropy.table.Table`
+
+        """
+        kwargs = {}
+        kwargs["creation"] = CreatorMetaData()
+        telescope = np.atleast_1d(table["TELESCOP"].data[0])
+        obs_id = np.atleast_1d(table["OBS_ID"].data[0].astype(str))
+        observation_mode = np.atleast_1d(table["OBS_MODE"].data[0])
+
+        obs_info = []
+        for i in range(len(obs_id)):
+            obs_meta = ObsInfoMetaData(
+                **{
+                    "telescope": telescope[i],
+                    "obs_id": obs_id[i],
+                    "observation_mode": observation_mode[i],
+                }
+            )
+            obs_info.append(obs_meta)
+        kwargs["obs_info"] = obs_info
+
+        pointing_radec, pointing_altaz = None, None
+        if "RA_PNT" in table.colnames:
+            pointing_radec = SkyCoord(
+                ra=table["RA_PNT"].data[0], dec=table["DEC_PNT"].data[0], unit="deg"
+            )
+        if "ALT_PNT" in table.colnames:
+            pointing_altaz = SkyCoord(
+                alt=table["ALT_PNT"].data[0],
+                az=table["AZ_PNT"].data[0],
+                unit="deg",
+                frame="altaz",
+            )
+        pointings = []
+        for pra, paz in zip(
+            np.atleast_1d(pointing_radec), np.atleast_1d(pointing_altaz)
+        ):
+            pointings.append(PointingInfoMetaData(radec_mean=pra, altaz_mean=paz))
+        kwargs["pointing"] = pointings
+        return cls(**kwargs)

gammapy/datasets/tests/test_map.py

@@ -602,6 +602,7 @@ def test_map_dataset_fits_io(tmp_path, sky_model, geom, geom_etrue):
     dataset_new = MapDataset.read(tmp_path / "test.fits")
 
     assert dataset_new.name == "test"
+    assert_allclose(dataset.meta.creation.date.mjd, dataset_new.meta.creation.date.mjd)
 
     assert dataset_new.mask.data.dtype == bool
 

gammapy/datasets/tests/test_metadata.py

@@ -0,0 +1,111 @@
+import pytest
+from numpy.testing import assert_allclose
+from astropy.coordinates import SkyCoord
+from pydantic import ValidationError
+from gammapy.datasets import MapDatasetMetaData
+from gammapy.utils.metadata import ObsInfoMetaData, PointingInfoMetaData
+
+
+def test_meta_default():
+    meta = MapDatasetMetaData()
+    assert meta.creation.creator.split()[0] == "Gammapy"
+    assert meta.obs_info is None
+
+
+def test_mapdataset_metadata():
+    position = SkyCoord(83.6287, 22.0147, unit="deg", frame="icrs")
+    obs_info_input = {
+        "telescope": "cta-north",
+        "instrument": "lst",
+        "observation_mode": "wobble",
+        "obs_id": 112,
+    }
+    input = {
+        "obs_info": ObsInfoMetaData(**obs_info_input),
+        "pointing": PointingInfoMetaData(radec_mean=position),
+        "optional": dict(test=0.5, other=True),
+    }
+    meta = MapDatasetMetaData(**input)
+
+    assert meta.obs_info.telescope == "cta-north"
+    assert meta.obs_info.instrument == "lst"
+    assert meta.obs_info.observation_mode == "wobble"
+    assert_allclose(meta.pointing.radec_mean.dec.value, 22.0147)
+    assert_allclose(meta.pointing.radec_mean.ra.deg, 83.6287)
+    assert meta.obs_info.obs_id == 112
+    assert meta.optional["other"] is True
+    assert meta.creation.creator.split()[0] == "Gammapy"
+    assert meta.event_type is None
+
+    with pytest.raises(ValidationError):
+        meta.pointing = 2.0
+
+    input_bad = input.copy()
+    input_bad["bad"] = position
+
+    with pytest.raises(ValueError):
+        MapDatasetMetaData(**input_bad)
+
+
+def test_mapdataset_metadata_lists():
+    obs_info_input1 = {
+        "telescope": "cta-north",
+        "instrument": "lst",
+        "observation_mode": "wobble",
+        "obs_id": 111,
+    }
+    obs_info_input2 = {
+        "telescope": "cta-north",
+        "instrument": "lst",
+        "observation_mode": "wobble",
+        "obs_id": 112,
+    }
+    input = {
+        "obs_info": [
+            ObsInfoMetaData(**obs_info_input1),
+            ObsInfoMetaData(**obs_info_input2),
+        ],
+        "pointing": [
+            PointingInfoMetaData(
+                radec_mean=SkyCoord(83.6287, 22.0147, unit="deg", frame="icrs")
+            ),
+            PointingInfoMetaData(
+                radec_mean=SkyCoord(83.1287, 22.5147, unit="deg", frame="icrs")
+            ),
+        ],
+    }
+    meta = MapDatasetMetaData(**input)
+    assert meta.obs_info[0].telescope == "cta-north"
+    assert meta.obs_info[0].instrument == "lst"
+    assert meta.obs_info[0].observation_mode == "wobble"
+    assert_allclose(meta.pointing[0].radec_mean.dec.value, 22.0147)
+    assert_allclose(meta.pointing[1].radec_mean.ra.deg, 83.1287)
+    assert meta.obs_info[0].obs_id == 111
+    assert meta.obs_info[1].obs_id == 112
+    assert meta.optional is None
+    assert meta.event_type is None
+
+
+def test_mapdataset_metadata_stack():
+    input1 = {
+        "obs_info": ObsInfoMetaData(**{"obs_id": 111}),
+        "pointing": PointingInfoMetaData(
+            radec_mean=SkyCoord(83.6287, 22.5147, unit="deg", frame="icrs")
+        ),
+        "optional": dict(test=0.5, other=True),
+    }
+
+    input2 = {
+        "obs_info": ObsInfoMetaData(**{"instrument": "H.E.S.S.", "obs_id": 112}),
+        "pointing": PointingInfoMetaData(
+            radec_mean=SkyCoord(83.6287, 22.0147, unit="deg", frame="icrs")
+        ),
+        "optional": dict(test=0.1, other=False),
+    }
+
+    meta1 = MapDatasetMetaData(**input1)
+    meta2 = MapDatasetMetaData(**input2)
+
+    meta = meta1.stack(meta2)
+    assert meta.creation.creator.split()[0] == "Gammapy"
+    assert meta.obs_info is None

gammapy/makers/map.py

@@ -3,6 +3,7 @@
 import astropy.units as u
 from astropy.table import Table
 from regions import PointSkyRegion
+from gammapy.datasets import MapDatasetMetaData
 from gammapy.irf import EDispKernelMap, PSFMap, RecoPSFMap
 from gammapy.maps import Map
 from .core import Maker
@@ -361,6 +362,10 @@ def make_meta_table(observation):
 
         return meta_table
 
+    @staticmethod
+    def _make_metadata(table):
+        return MapDatasetMetaData._from_meta_table(table)
+
     def run(self, dataset, observation):
         """Make map dataset.
 
@@ -378,6 +383,7 @@ def run(self, dataset, observation):
         """
         kwargs = {"gti": observation.gti}
         kwargs["meta_table"] = self.make_meta_table(observation)
+        kwargs["meta"] = self._make_metadata(kwargs["meta_table"])
 
         mask_safe = Map.from_geom(dataset.counts.geom, dtype=bool)
         mask_safe.data[...] = True