Merge pull request #62 from nbiederbeck/skymap-dl3

Skymap dl3

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))

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}")
 

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))

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"]
 ]
 
 

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}"

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}"