diff --git a/scripts/calc_skymap_gammas.py b/scripts/calc_skymap_gammas.py
new file mode 100644
index 0000000..34a1257
--- /dev/null
+++ b/scripts/calc_skymap_gammas.py
@@ -0,0 +1,76 @@
+import json
+from argparse import ArgumentParser
+
+import numpy as np
+from astropy import units as u
+from astropy.coordinates import SkyCoord
+from gammapy.data import DataStore
+from gammapy.maps import MapAxis, WcsGeom, WcsNDMap
+
+parser = ArgumentParser()
+parser.add_argument("-i", "--input-path", required=True)
+parser.add_argument("-o", "--output-path", required=True)
+parser.add_argument("-c", "--config", required=True)
+parser.add_argument("--obs-id", required=True)
+parser.add_argument(
+    "--width",
+    help="Width of skymap",
+    default="3 deg",
+    type=u.Quantity,
+)
+parser.add_argument("--n-bins", default=100)
+args = parser.parse_args()
+
+
+@u.quantity_input(width=u.deg)
+def main(input_path, config, output_path, obs_id, width, n_bins):  # noqa: PLR0913
+    with open(config, "r") as f:
+        config = json.load(f)
+
+    ds = DataStore.from_dir(input_path)
+
+    obs = ds.obs(int(obs_id), ["aeff"])
+
+    events = obs.events
+
+    evts = events.radec
+
+    source = SkyCoord(
+        config["DataReductionFITSWriter"]["source_ra"],
+        config["DataReductionFITSWriter"]["source_dec"],
+        frame="icrs",
+    )
+
+    r = width / 2
+    bins = u.Quantity(
+        [
+            np.linspace(source.ra - r, source.ra + r, n_bins + 1),
+            np.linspace(source.dec - r, source.dec + r, n_bins + 1),
+        ],
+    )
+
+    hist, *edges = np.histogram2d(
+        evts.ra.to_value(u.deg),
+        evts.dec.to_value(u.deg),
+        bins.to_value(u.deg),
+    )
+
+    energy_axis = MapAxis.from_edges(
+        u.Quantity([1e-5, 1e5], u.TeV),  # basically [-inf, inf], but finite
+        name="energy",
+    )
+    geom = WcsGeom.create(
+        (n_bins, n_bins),
+        skydir=source,
+        binsz=width / n_bins,
+        width=width,
+        axes=[energy_axis],
+    )
+
+    skymap = WcsNDMap(geom, hist)
+
+    skymap.write(output_path, overwrite=True)
+
+
+if __name__ == "__main__":
+    main(**vars(args))
diff --git a/scripts/plot_skymap.py b/scripts/plot_skymap.py
index 8e10641..a16b371 100644
--- a/scripts/plot_skymap.py
+++ b/scripts/plot_skymap.py
@@ -27,7 +27,7 @@ def main(input_path, output_path):
 
     fig, ax = plt.subplots()
 
-    ax.pcolormesh(
+    mesh = ax.pcolormesh(
         *edges.to_value(angle),
         skymap.data[0, ...],
         rasterized=True,
@@ -41,6 +41,8 @@ def main(input_path, output_path):
         label="Source",
     )
 
+    fig.colorbar(mesh, ax=ax)
+
     ax.set_xlabel(f"RA / {angle}")
     ax.set_ylabel(f"Dec / {angle}")
 
diff --git a/scripts/plot_skymap_dl3.py b/scripts/plot_skymap_dl3.py
new file mode 100644
index 0000000..5a616e1
--- /dev/null
+++ b/scripts/plot_skymap_dl3.py
@@ -0,0 +1,95 @@
+from argparse import ArgumentParser
+
+import numpy as np
+from astropy import units as u
+from gammapy.maps import WcsNDMap
+from matplotlib import pyplot as plt
+from matplotlib.collections import PatchCollection
+from matplotlib.patches import Circle
+
+parser = ArgumentParser()
+parser.add_argument("-i", "--input-path", required=True)
+parser.add_argument("-o", "--output-path", required=True)
+parser.add_argument(
+    "--width",
+    help="Width of skymap",
+    default="4 deg",
+    type=u.Quantity,
+)
+parser.add_argument("--n-bins", default=100)
+args = parser.parse_args()
+
+
+@u.quantity_input(width=u.deg)
+def main(input_path, output_path, width, n_bins):
+    skymap = WcsNDMap.read(input_path)
+
+    geom = skymap.geom
+    source = geom.center_skydir
+
+    edges = u.Quantity(
+        [
+            u.Quantity(geom.pix_to_coord([i, j]))
+            for i, j in zip(range(int(geom.npix[0])), range(int(geom.npix[1])))
+        ],
+    ).T
+
+    patches = PatchCollection(
+        [
+            Circle(
+                (source.ra.to_value(u.deg), source.dec.to_value(u.deg)),
+                radius=0.3,
+                color="k",
+                ec="k",
+                fc="#0000",
+                fill=False,
+            ),
+        ]
+        + [
+            Circle(
+                (
+                    source.ra.to_value(u.deg) + 0.8 * np.sin(phi),
+                    source.dec.to_value(u.deg) + 0.8 * np.cos(phi),
+                ),
+                radius=0.3,
+                color="w",
+                ec="w",
+                fc="#0000",
+                fill=False,
+            )
+            for phi in np.linspace(0, 2 * np.pi, 7)
+        ],
+        match_original=True,
+    )
+
+    fig, ax = plt.subplots()
+
+    mesh = ax.pcolormesh(
+        *edges.to_value(u.deg),
+        skymap.data[0, ...],
+        rasterized=True,
+    )
+
+    fig.colorbar(mesh, ax=ax)
+
+    ax.set_xlabel("RA / deg")
+    ax.set_ylabel("Dec / deg")
+    ax.set_aspect(1)
+
+    ax.scatter(
+        source.ra.to_value(u.deg),
+        source.dec.to_value(u.deg),
+        color="w",
+        edgecolor="k",
+        label="Source",
+    )
+
+    ax.add_collection(patches)
+
+    ax.legend()
+
+    fig.savefig(output_path)
+
+
+if __name__ == "__main__":
+    main(**vars(args))
diff --git a/workflow/Snakefile b/workflow/Snakefile
index 6a250a9..8342ea5 100644
--- a/workflow/Snakefile
+++ b/workflow/Snakefile
@@ -29,7 +29,7 @@ analyses = [
 dl3_plots = [
     build_dir / f"plots/{plot}/{run_id}.pdf"
     for run_id in RUN_IDS + ["stacked"]
-    for plot in ["theta2", "skymap"]
+    for plot in ["theta2", "skymap", "skymap-dl3"]
 ]
 
 
diff --git a/workflow/rules/create_dl3.smk b/workflow/rules/create_dl3.smk
index 6c93972..638f1f1 100644
--- a/workflow/rules/create_dl3.smk
+++ b/workflow/rules/create_dl3.smk
@@ -179,3 +179,18 @@ rule stack_skymaps:
         script="scripts/stack_skymap.py",
     shell:
         "python {input.script} -i {input.data} -o {output}"
+
+
+rule stack_skymaps_dl3:
+    conda:
+        lstchain_env
+    output:
+        build_dir / "dl3/skymap-dl3/stacked.fits",
+    input:
+        data=expand(
+            build_dir / "dl3/skymap-dl3/{run_id}.fits",
+            run_id=RUN_IDS,
+        ),
+        script="scripts/stack_skymap.py",
+    shell:
+        "python {input.script} -i {input.data} -o {output}"
diff --git a/workflow/rules/post_dl3.smk b/workflow/rules/post_dl3.smk
index bc005ee..b4416b0 100644
--- a/workflow/rules/post_dl3.smk
+++ b/workflow/rules/post_dl3.smk
@@ -215,3 +215,37 @@ rule model_best_fit:
             --model-config {input.model} \
             -o {output} \
         """
+
+
+rule calc_skymap_per_obs:
+    output:
+        build_dir / "dl3/skymap-dl3/{run_id}.fits",
+    input:
+        data=build_dir / "dl3/dl3_LST-1.Run{run_id}.fits.gz",
+        script="scripts/calc_skymap_gammas.py",
+        config=irf_config_path,
+        index=build_dir / "dl3/hdu-index.fits.gz",
+    wildcard_constraints:
+        run_id="\d+",  # dont match on "stacked".
+    resources:
+        # mem_mb=16000,
+        time=5,
+    conda:
+        gammapy_env
+    shell:
+        "python {input.script} -i {build_dir}/dl3 -o {output} --obs-id {wildcards.run_id} --config {input.config}"
+
+
+rule plot_skymap_dl3:
+    output:
+        build_dir / "plots/skymap-dl3/{runid}.pdf",
+    input:
+        data=build_dir / "dl3/skymap-dl3/{runid}.fits",
+        script="scripts/plot_skymap_dl3.py",
+        rc=os.environ.get("MATPLOTLIBRC", config_dir / "matplotlibrc"),
+    conda:
+        gammapy_env
+    resources:
+        time=5,
+    shell:
+        "MATPLOTLIBRC={input.rc} python {input.script} -i {input.data} -o {output}"