Repository: gammapy/gammapy
Pull request #5875 (2025-05-19 21:37:07)

Export creation metadata on DL3 and DL4 products
<!-- These are hidden comments. You can leave them or delete them. -->

<!-- Don't be shy! Github issues or pull requests are welcome any time. -->
<!-- But if you first want to chat, use the Gammapy Slack. -->
<!-- See https://gammapy.org/contact.html -->

**Description**
<!-- What is the motivation for this pull request? -->
<!-- Briefly: what changes are done here? -->
<!-- If this is to address a Github issue, mention its number and use the specific keywords to create a link -->
<!-- See docs.github.com/en/issues/tracking-your-work-with-issues/linking-a-pull-request-to-an-issue -->
<!-- Example: This PR is to resolve #42 -->

This pull request introduces support for adding creation metadata on DL3 and DL4 objects serialized to FITS.

The proposed approach is to add a `creation` argument on `to_hdulist` functions of objects which do not contain regular metadata (such as `Map`) and to make sure the `CreatorMetaData` is correctly instantiated in all objects containing proper `MetaData` objects.    

This approach should make it possible for a user or gammapy dependent code to set/adapt creation information.
We could also introduce a global variable containing the default string for the `CREATOR` information.

This is a first step to solve #5598. The additional step will be to check DL5 products. 

**Dear reviewer**
<!-- Let the reviewer and Gammapy team know what you want: -->
<!-- * Is this ready for review? Or is it work in progress and you want some feedback? -->
<!-- * Do you want to go through review here? Or if someone should just finish this up and merge it? -->
<!-- Do you have any specific questions, e.g. about API or implementation? -->
<!-- Do you include a test executing the new code you're adding (to make sure it runs)? -->
<!-- Do you include some documentation? Is it needed? -->


registerrier: Decision taken during dev-call: do not introduce `creation` keyword into `Map.to_hdulist()`. This is too specific a keyword and adding the information on the maps HDUs can be done directly with the `MapDataset.to_hdulist()` method.

This can be revisited once we have defined a proper metadata scheme for `Map` objects.

Patch:
diff --git a/gammapy/data/__init__.py b/gammapy/data/__init__.py
index 28d91b1f8..a167f564c 100644
--- a/gammapy/data/__init__.py
+++ b/gammapy/data/__init__.py
@@ -1,12 +1,13 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """Data and observation handling."""
+
 from gammapy.utils.observers import observatory_locations
 from .data_store import DataStore
 from .event_list import EventList
 from .filters import ObservationFilter
 from .gti import GTI
 from .hdu_index_table import HDUIndexTable
-from .metadata import EventListMetaData, ObservationMetaData
+from .metadata import EventListMetaData, ObservationMetaData, GTIMetaData
 from .obs_table import ObservationTable
 from .observations import Observation, Observations
 from .pointing import FixedPointingInfo, PointingInfo, PointingMode
@@ -20,6 +21,7 @@ __all__ = [
     "ObservationMetaData",
     "FixedPointingInfo",
     "GTI",
+    "GTIMetaData",
     "HDUIndexTable",
     "Observation",
     "ObservationFilter",
diff --git a/gammapy/data/event_list.py b/gammapy/data/event_list.py
index 3dede50d6..cba70cf96 100644
--- a/gammapy/data/event_list.py
+++ b/gammapy/data/event_list.py
@@ -20,6 +20,7 @@ from gammapy.utils.testing import Checker
 from gammapy.utils.time import time_ref_from_dict
 from .metadata import EventListMetaData
 from gammapy.utils.deprecation import deprecated_renamed_argument
+from gammapy.utils.metadata import CreatorMetaData
 
 __all__ = ["EventList"]
 
@@ -92,7 +93,7 @@ class EventList:
 
     def __init__(self, table, meta=None):
         self.table = table
-        self.meta = meta
+        self.meta = meta or EventListMetaData()
 
     def _repr_html_(self):
         try:
@@ -148,7 +149,17 @@ class EventList:
         if format != "gadf":
             raise ValueError(f"Only the 'gadf' format supported, got {format}")
 
-        return fits.BinTableHDU(self.table, name="EVENTS")
+        bin_table = fits.BinTableHDU(self.table, name="EVENTS")
+
+        # A priori don't change creator information
+        if self.meta.creation is None:
+            self.meta.creation = CreatorMetaData()
+        else:
+            self.meta.creation.update_time()
+
+        bin_table.header.update(self.meta.to_header())
+
+        return bin_table
 
     # TODO: Pass metadata here. Also check that specific meta contents are consistent
     @classmethod
diff --git a/gammapy/data/gti.py b/gammapy/data/gti.py
index b6e63d932..c71dab820 100644
--- a/gammapy/data/gti.py
+++ b/gammapy/data/gti.py
@@ -8,8 +8,9 @@ import astropy.units as u
 from astropy.io import fits
 from astropy.table import Table, vstack
 from astropy.time import Time
+from gammapy.utils.metadata import CreatorMetaData
 from gammapy.utils.scripts import make_path
-from gammapy.utils.time import TIME_REF_DEFAULT, time_ref_from_dict, time_ref_to_dict
+from gammapy.utils.time import TIME_REF_DEFAULT, time_ref_from_dict
 from .metadata import GTIMetaData
 
 __all__ = ["GTI"]
@@ -201,10 +202,13 @@ class GTI:
             raise ValueError(f'Only the "gadf" format supported, got {format}')
 
         # Don't impose the scale. GADF does not require it to be TT
-        meta = time_ref_to_dict(self.time_ref, scale=self.time_ref.scale)
+        if self._meta.creation is None:
+            self._meta.creation = CreatorMetaData()
+        hdr = self._meta.to_header(format)
+
         start = self.time_start - self.time_ref
         stop = self.time_stop - self.time_ref
-        table = Table({"START": start.to("s"), "STOP": stop.to("s")}, meta=meta)
+        table = Table({"START": start.to("s"), "STOP": stop.to("s")}, meta=hdr)
 
         return fits.BinTableHDU(table, name="GTI")
 
diff --git a/gammapy/data/metadata.py b/gammapy/data/metadata.py
index c5033e439..693c3b3e0 100644
--- a/gammapy/data/metadata.py
+++ b/gammapy/data/metadata.py
@@ -11,6 +11,7 @@ from gammapy.utils.metadata import (
     TargetMetaData,
     TimeInfoMetaData,
 )
+from gammapy.utils.time import time_ref_from_dict, time_ref_to_dict
 from gammapy.utils.types import EarthLocationType, TimeType
 
 __all__ = ["ObservationMetaData", "GTIMetaData", "EventListMetaData"]
@@ -38,6 +39,7 @@ EVENTLIST_METADATA_FITS_KEYS = {
 }
 
 METADATA_FITS_KEYS["eventlist"] = EVENTLIST_METADATA_FITS_KEYS
+METADATA_FITS_KEYS["GTI"] = {}
 
 
 class ObservationMetaData(MetaData):
@@ -126,13 +128,29 @@ class GTIMetaData(MetaData):
     ----------
     reference_time : Time, str
         The GTI reference time.
+    creation : `~gammapy.utils.metadata.CreatorMetaData`
+        The creation metadata.
     """
 
     _tag: ClassVar[Literal["GTI"]] = "GTI"
     reference_time: Optional[TimeType] = None
+    creation: Optional[CreatorMetaData] = None
 
+    def to_header(self, format="gadf"):
+        result = super().to_header(format)
+        if self.reference_time is not None:
+            result.update(
+                time_ref_to_dict(self.reference_time, scale=self.reference_time.scale)
+            )
+        return result
+
+    @classmethod
     def from_header(cls, header, format="gadf"):
         meta = super().from_header(header, format)
+        try:
+            meta.reference_time = time_ref_from_dict(header)
+        except KeyError:
+            meta.reference_time = None
 
         return meta
 
diff --git a/gammapy/datasets/io.py b/gammapy/datasets/io.py
index d669a76c4..d57bdb1fb 100644
--- a/gammapy/datasets/io.py
+++ b/gammapy/datasets/io.py
@@ -10,6 +10,7 @@ from gammapy.irf import EDispKernel, EDispKernelMap, PSFMap
 from gammapy.maps import RegionNDMap, Map
 from gammapy.modeling.models import create_fermi_isotropic_diffuse_model, Models
 from gammapy.utils.scripts import read_yaml, make_path
+from gammapy.utils.metadata import CreatorMetaData
 from .spectrum import SpectrumDatasetOnOff
 from .utils import create_map_dataset_from_dl4
 
@@ -73,11 +74,15 @@ class OGIPDatasetWriter(DatasetWriter):
     checksum : bool
         When True adds both DATASUM and CHECKSUM cards to the headers written to the files.
         Default is False.
+    creation : `~gammapy.utils.metadata.CreatorMetaData`, optional.
+        The creation metadata to write to disk. If None, use default creator metadata object.
     """
 
     tag = ["ogip", "ogip-sherpa"]
 
-    def __init__(self, filename, format="ogip", overwrite=False, checksum=False):
+    def __init__(
+        self, filename, format="ogip", overwrite=False, checksum=False, creation=None
+    ):
         filename = make_path(filename)
         filename.parent.mkdir(exist_ok=True, parents=True)
 
@@ -85,6 +90,7 @@ class OGIPDatasetWriter(DatasetWriter):
         self.format = format
         self.overwrite = overwrite
         self.checksum = checksum
+        self.creation = creation or CreatorMetaData()
 
     @staticmethod
     def get_filenames(filename):
@@ -150,6 +156,8 @@ class OGIPDatasetWriter(DatasetWriter):
         """
         filenames = self.get_filenames(self.filename)
 
+        self.creation.update_time()
+
         self.write_pha(dataset, filename=self.filename)
 
         path = self.filename.parent
@@ -191,12 +199,14 @@ class OGIPDatasetWriter(DatasetWriter):
 
         """
         aeff = dataset.exposure / dataset.exposure.meta["livetime"]
-        aeff.write(
-            filename=filename,
-            overwrite=self.overwrite,
-            format=self.format.replace("ogip", "ogip-arf"),
-            checksum=self.checksum,
-        )
+
+        filename = make_path(filename)
+        hdulist = aeff.to_hdulist(format=self.format.replace("ogip", "ogip-arf"))
+
+        for hdu in hdulist:
+            hdu.header.update(self.creation.to_header())
+
+        hdulist.writeto(filename, overwrite=self.overwrite, checksum=self.checksum)
 
     def to_counts_hdulist(self, dataset, is_bkg=False):
         """Convert counts region map to hdulist.
@@ -230,6 +240,9 @@ class OGIPDatasetWriter(DatasetWriter):
         # adapt meta data
         table.meta.update(meta)
         hdulist["SPECTRUM"] = fits.BinTableHDU(table)
+
+        for hdu in hdulist:
+            hdu.header.update(self.creation.to_header())
         return hdulist
 
     def write_pha(self, dataset, filename):
diff --git a/gammapy/datasets/map.py b/gammapy/datasets/map.py
index a0cce86a6..1fc11c7e4 100644
--- a/gammapy/datasets/map.py
+++ b/gammapy/datasets/map.py
@@ -1518,6 +1518,8 @@ class MapDataset(Dataset):
         header = hdu_primary.header
         header["NAME"] = self.name
         header.update(self.meta.to_header())
+        creation = self.meta.creation
+        creation.update_time()
 
         hdulist = fits.HDUList([hdu_primary])
         if self.counts is not None:
@@ -1547,6 +1549,9 @@ class MapDataset(Dataset):
         if self.meta_table is not None:
             hdulist.append(fits.BinTableHDU(self.meta_table, name="META_TABLE"))
 
+        for hdu in hdulist:
+            hdu.header.update(creation.to_header())
+
         return hdulist
 
     @classmethod
@@ -2983,6 +2988,8 @@ class MapDatasetOnOff(MapDataset):
         hdulist = super().to_hdulist()
         exclude_primary = slice(1, None)
 
+        creation = self.meta.creation
+
         del hdulist["BACKGROUND"]
         del hdulist["BACKGROUND_BANDS"]
 
@@ -2997,6 +3004,9 @@ class MapDatasetOnOff(MapDataset):
                 exclude_primary
             ]
 
+        for hdu in hdulist:
+            hdu.header.update(creation.to_header())
+
         return hdulist
 
     @classmethod
diff --git a/gammapy/datasets/spectrum.py b/gammapy/datasets/spectrum.py
index a4f218d5e..fdf1fd7e2 100644
--- a/gammapy/datasets/spectrum.py
+++ b/gammapy/datasets/spectrum.py
@@ -3,6 +3,7 @@ import logging
 import matplotlib.pyplot as plt
 from matplotlib.gridspec import GridSpec
 from gammapy.utils.scripts import make_path
+from gammapy.utils.metadata import CreatorMetaData
 from .map import MapDataset, MapDatasetOnOff
 from .utils import get_axes
 
@@ -357,8 +358,13 @@ class SpectrumDatasetOnOff(PlotMixin, MapDatasetOnOff):
         if format == "gadf":
             super().write(filename=filename, overwrite=overwrite, checksum=checksum)
         elif format in ["ogip", "ogip-sherpa"]:
+            creation = self.meta.creation or CreatorMetaData()
             writer = OGIPDatasetWriter(
-                filename=filename, format=format, overwrite=overwrite, checksum=checksum
+                filename=filename,
+                format=format,
+                overwrite=overwrite,
+                checksum=checksum,
+                creation=creation,
             )
             writer.write(self)
         else:
diff --git a/gammapy/irf/edisp/kernel.py b/gammapy/irf/edisp/kernel.py
index 26798cd78..55abcfeb5 100644
--- a/gammapy/irf/edisp/kernel.py
+++ b/gammapy/irf/edisp/kernel.py
@@ -10,6 +10,7 @@ from matplotlib.colors import PowerNorm
 from gammapy.maps import MapAxis
 from gammapy.maps.axes import UNIT_STRING_FORMAT
 from gammapy.utils.scripts import make_path
+from gammapy.utils.metadata import CreatorMetaData
 from gammapy.visualization.utils import add_colorbar
 from ..core import IRF
 
@@ -310,13 +311,16 @@ class EDispKernel(IRF):
         else:
             raise ValueError(f"Unrecognized format: {format}")
 
-    def to_hdulist(self, format="ogip", **kwargs):
+    def to_hdulist(self, format="ogip", creation=None, **kwargs):
         """Convert RMF to FITS HDU list format.
 
         Parameters
         ----------
         format : {"ogip", "ogip-sherpa"}
             Format to use. Default is "ogip".
+        creation : `~gammapy.utils.metadata.CreatorMetadata`, optional
+            Creation metadata to add to the file. If None, default metadata is added.
+            Default is None.
 
         Returns
         -------
@@ -357,6 +361,12 @@ class EDispKernel(IRF):
         ebounds_hdu = self.axes["energy"].to_table_hdu(format=format)
         prim_hdu = fits.PrimaryHDU()
 
+        creation = creation or CreatorMetaData()
+        creation.update_time()
+
+        for hd in [prim_hdu, hdu, ebounds_hdu]:
+            hd.header.update(creation.to_header())
+
         return fits.HDUList([prim_hdu, hdu, ebounds_hdu])
 
     def to_table(self, format="ogip"):
@@ -429,7 +439,7 @@ class EDispKernel(IRF):
 
         return table
 
-    def write(self, filename, format="ogip", checksum=False, **kwargs):
+    def write(self, filename, format="ogip", checksum=False, creation=None, **kwargs):
         """Write to file.
 
         Parameters
@@ -440,10 +450,13 @@ class EDispKernel(IRF):
             Format to use. Default is "ogip".
         checksum : bool
             If True checks both DATASUM and CHECKSUM cards in the file headers. Default is False.
+        creation : `~gammapy.utils.metadata.CreatorMetadata`, optional
+            Creation metadata to add to the file. If None, default metadata is added.
+            Default is None.
 
         """
         filename = str(make_path(filename))
-        hdulist = self.to_hdulist(format=format)
+        hdulist = self.to_hdulist(format=format, creation=creation)
 
         hdulist.writeto(filename, checksum=checksum, **kwargs)
 
diff --git a/gammapy/maps/region/ndmap.py b/gammapy/maps/region/ndmap.py
index d3497b0e6..32cc73333 100644
--- a/gammapy/maps/region/ndmap.py
+++ b/gammapy/maps/region/ndmap.py
@@ -513,7 +513,12 @@ class RegionNDMap(Map):
             )
 
     def write(
-        self, filename, overwrite=False, format="gadf", hdu="SKYMAP", checksum=False
+        self,
+        filename,
+        overwrite=False,
+        format="gadf",
+        hdu="SKYMAP",
+        checksum=False,
     ):
         """Write map to file.
 
diff --git a/gammapy/utils/metadata.py b/gammapy/utils/metadata.py
index 3409604c0..d34422451 100644
--- a/gammapy/utils/metadata.py
+++ b/gammapy/utils/metadata.py
@@ -196,6 +196,10 @@ class CreatorMetaData(MetaData):
     date: Optional[TimeType] = Field(default_factory=Time.now)
     origin: Optional[str] = None
 
+    def update_time(self):
+        """Change creation date to Time.now()"""
+        self.date = Time.now()
+
 
 class ObsInfoMetaData(MetaData):
     """General metadata information about the observation.


Test patch:
diff --git a/gammapy/data/tests/test_event_list.py b/gammapy/data/tests/test_event_list.py
index f6d714f21..e6cf6a6e8 100644
--- a/gammapy/data/tests/test_event_list.py
+++ b/gammapy/data/tests/test_event_list.py
@@ -80,6 +80,12 @@ class TestEventListBase:
             obs = Observation(events=self.events, gti=gti.table, pointing=pointing)
             obs.write("test.fits", overwrite=True)
 
+    def test_eventlist_hdu_creation_metadata(self):
+        hdu = self.events.to_table_hdu(format="gadf")
+        assert "CREATOR" in hdu.header
+        assert "CREATED" in hdu.header
+        assert hdu.header["CREATOR"] == "SASH FITS::EventListWriter"
+
 
 @requires_data()
 class TestEventListHESS:
diff --git a/gammapy/data/tests/test_gti.py b/gammapy/data/tests/test_gti.py
index 8bf6e2b08..164dd4afc 100644
--- a/gammapy/data/tests/test_gti.py
+++ b/gammapy/data/tests/test_gti.py
@@ -264,3 +264,11 @@ def test_gti_from_time():
 
     assert_time_allclose(gti.table["START"], start)
     assert_time_allclose(gti.table["STOP"], stop)
+
+
+def test_gti_creation_metadata():
+    time_ref = Time("2010-01-01")
+    gti = make_gti({"START": [0, 2] * u.s, "STOP": [1, 3] * u.s}, time_ref=time_ref)
+
+    hdu = gti.to_table_hdu(format="gadf")
+    assert "MJDREFI" in hdu.header
diff --git a/gammapy/data/tests/test_metadata.py b/gammapy/data/tests/test_metadata.py
index b0c6b3dd2..2d3660365 100644
--- a/gammapy/data/tests/test_metadata.py
+++ b/gammapy/data/tests/test_metadata.py
@@ -4,8 +4,13 @@ from numpy.testing import assert_allclose
 from astropy.coordinates import SkyCoord
 from astropy.io import fits
 from pydantic import ValidationError
-from gammapy.data import EventListMetaData, ObservationMetaData
-from gammapy.utils.metadata import ObsInfoMetaData, PointingInfoMetaData, TargetMetaData
+from gammapy.data import EventListMetaData, ObservationMetaData, GTIMetaData
+from gammapy.utils.metadata import (
+    ObsInfoMetaData,
+    PointingInfoMetaData,
+    TargetMetaData,
+    CreatorMetaData,
+)
 from gammapy.utils.scripts import make_path
 from gammapy.utils.testing import requires_data
 
@@ -114,3 +119,23 @@ def test_eventlist_metadata():
     input_bad["location"] = "bad"
     with pytest.raises(ValueError):
         EventListMetaData(**input_bad)
+
+
+def test_gti_metadata():
+    input = {
+        "reference_time": "2023-01-01 00:00:00",
+        "creation": CreatorMetaData(creator="Test"),
+    }
+    meta = GTIMetaData(**input)
+
+    hdr = meta.to_header()
+
+    assert meta.creation.creator == "Test"
+    assert_allclose(meta.reference_time.mjd, 59945.0)
+
+    assert hdr["MJDREFI"] == 59945
+    assert hdr["TIMESYS"] == "utc"
+    assert hdr["CREATOR"] == "Test"
+
+    meta_new = GTIMetaData.from_header(hdr)
+    assert_allclose(meta_new.reference_time.mjd, 59945.0)
diff --git a/gammapy/datasets/tests/test_map.py b/gammapy/datasets/tests/test_map.py
index 3b0e6aeba..d9bd97069 100644
--- a/gammapy/datasets/tests/test_map.py
+++ b/gammapy/datasets/tests/test_map.py
@@ -8,6 +8,7 @@ import astropy.units as u
 from astropy.coordinates import SkyCoord
 from astropy.io import fits
 from astropy.table import Table
+from astropy.time import Time
 from astropy.utils.exceptions import AstropyUserWarning
 from regions import CircleSkyRegion
 import gammapy.irf.psf.map as psf_map_module
@@ -746,6 +747,33 @@ def test_map_dataset_fits_io(tmp_path, sky_model, geom, geom_etrue):
     assert_allclose(stacked1.edisp.edisp_map, stacked.edisp.edisp_map)
 
 
+@requires_data()
+def test_map_dataset_fits_creation_metadata(tmp_path, sky_model, geom, geom_etrue):
+    dataset = get_map_dataset(geom, geom_etrue)
+
+    bkg_model = FoVBackgroundModel(dataset_name=dataset.name)
+    dataset.models = [sky_model, bkg_model]
+
+    dataset.meta.creation.creator = "MySoftware"
+    dataset.meta.creation.origin = "MyOrganization"
+
+    dataset.counts = dataset.npred()
+    dataset.mask_safe = dataset.mask_fit
+    gti = GTI.create([0 * u.s], [1 * u.h], reference_time="2010-01-01T00:00:00")
+    dataset.gti = gti
+
+    dataset.write(tmp_path / "test.fits")
+
+    hdul = fits.open(tmp_path / "test.fits")
+    for hdu in hdul:
+        assert "CREATOR" in hdu.header
+        assert "CREATED" in hdu.header
+        assert hdu.header["CREATOR"] == "MySoftware"
+        assert hdu.header["ORIGIN"] == "MyOrganization"
+        # Check that day is OK
+        assert hdu.header["CREATED"][:10] == Time.now().iso[:10]
+
+
 @requires_data()
 def test_map_auto_psf_upsampling(sky_model, geom, geom_etrue):
     dataset_2 = get_map_dataset(geom, geom_etrue, name="test-2")
diff --git a/gammapy/datasets/tests/test_spectrum.py b/gammapy/datasets/tests/test_spectrum.py
index 0897dc6e7..f9fbdf21b 100644
--- a/gammapy/datasets/tests/test_spectrum.py
+++ b/gammapy/datasets/tests/test_spectrum.py
@@ -722,6 +722,16 @@ class TestSpectrumOnOff:
         assert_allclose(np.squeeze(grouped.acceptance), 4)
         assert_allclose(np.squeeze(grouped.acceptance_off), 40)
 
+    def test_spectrum_dataset_ogip_creation_metadata(self, tmp_path):
+        dataset = self.dataset.copy(name="test")
+        dataset.write(tmp_path / "test.fits", format="ogip")
+
+        for name in ["test.fits", "test_arf.fits", "test_bkg.fits", "test_rmf.fits"]:
+            with fits.open(tmp_path / name) as hdul:
+                for hdu in hdul[1:]:
+                    assert "CREATOR" in hdu.header
+                    assert "CREATED" in hdu.header
+
 
 @requires_data()
 class TestSpectralFit: