Add support for skymap with planned wobble pointings

docs/source/changes/29.added.rst

@@ -0,0 +1 @@
+Added support for plotting a skymap with the planned wobble pointings.

docs/source/example.ipynb

@@ -29,13 +29,14 @@
     "from iact_estimator import RESOURCES_PATH\n",
     "from iact_estimator.io import read_yaml\n",
     "from iact_estimator.core import (\n",
-    "\n",
     "    initialize_model,\n",
     "    prepare_data,\n",
     "    source_detection,\n",
     "    calculate,\n",
     ")\n",
-    "from iact_estimator.plots import plot_spectrum, plot_sed, plot_transit, plot_altitude_airmass"
+    "from iact_estimator.plots.physics import plot_spectrum, plot_sed\n",
+    "from iact_estimator.plots.observability import plot_transit, plot_altitude_airmass\n",
+    "from iact_estimator.plots.wobble_skymap import plot_skymap_with_wobbles, load_wobbles"
    ]
   },
   {
@@ -46,8 +47,8 @@
    },
    "outputs": [],
    "source": [
-    "output_path=Path.cwd()\n",
-    "config = read_yaml(RESOURCES_PATH/\"config.yml\")\n",
+    "output_path = Path.cwd()\n",
+    "config = read_yaml(RESOURCES_PATH / \"config.yml\")\n",
     "source_name = \"Crab\"\n",
     "\n",
     "observer = Observer.at_site(\"Roque de los Muchachos\")\n",
@@ -97,26 +98,26 @@
     "observer = Observer.at_site(\"Roque de los Muchachos\")\n",
     "\n",
     "date_time = DatetimePicker(\n",
-    "    value=datetime.now(timezone.utc),\n",
-    "    description='Select a datetime',\n",
-    "    disabled=False\n",
+    "    value=datetime.now(timezone.utc), description=\"Select a datetime\", disabled=False\n",
     ")\n",
     "\n",
     "crab = FixedTarget.from_name(\"Crab\")\n",
     "plot_crab = True if (crab.coord == target_source.coord) else False\n",
     "\n",
+    "\n",
     "def interactive_plot_transit(date_time):\n",
-    "        with quantity_support():\n",
-    "                plot_transit(\n",
-    "                        config,\n",
-    "                        source_name,\n",
-    "                        target_source,\n",
-    "                        observer,\n",
-    "                        time = Time(date_time).utc,\n",
-    "                        merge_profiles=True,\n",
-    "                        plot_crab=False,\n",
-    "                        savefig=False,\n",
-    "                )\n",
+    "    with quantity_support():\n",
+    "        plot_transit(\n",
+    "            config,\n",
+    "            source_name,\n",
+    "            target_source,\n",
+    "            observer,\n",
+    "            time=Time(date_time).utc,\n",
+    "            merge_profiles=True,\n",
+    "            plot_crab=False,\n",
+    "            savefig=False,\n",
+    "        )\n",
+    "\n",
     "\n",
     "interact(interactive_plot_transit, date_time=date_time)\n",
     "plt.show()"
@@ -136,25 +137,24 @@
    "outputs": [],
    "source": [
     "date_time = DatetimePicker(\n",
-    "    value=datetime.now(timezone.utc),\n",
-    "    description='Select a datetime',\n",
-    "    disabled=False\n",
+    "    value=datetime.now(timezone.utc), description=\"Select a datetime\", disabled=False\n",
     ")\n",
     "\n",
-    "def plot_alt(date_time):\n",
     "\n",
+    "def plot_alt(date_time):\n",
     "    print(date_time)\n",
     "\n",
     "    plot_altitude_airmass(\n",
-    "                config,\n",
-    "                source_name,\n",
-    "                target_source,\n",
-    "                observer,\n",
-    "                time=Time(date_time).utc,\n",
-    "                brightness_shading=True,\n",
-    "                airmass_yaxis=True,\n",
-    "                savefig=False,\n",
-    "            )\n",
+    "        config,\n",
+    "        source_name,\n",
+    "        target_source,\n",
+    "        observer,\n",
+    "        time=Time(date_time).utc,\n",
+    "        brightness_shading=True,\n",
+    "        airmass_yaxis=True,\n",
+    "        savefig=False,\n",
+    "    )\n",
+    "\n",
     "\n",
     "interact(plot_alt, date_time=date_time)\n",
     "plt.show()"
@@ -174,13 +174,13 @@
    "outputs": [],
    "source": [
     "plot_spectrum(\n",
-    "            config,\n",
-    "        plot_energy_bounds,\n",
-    "        assumed_spectrum,\n",
-    "        source_name,\n",
-    "        plotting_options,\n",
-    "        savefig=False,\n",
-    "        )"
+    "    config,\n",
+    "    plot_energy_bounds,\n",
+    "    assumed_spectrum,\n",
+    "    source_name,\n",
+    "    plotting_options,\n",
+    "    savefig=False,\n",
+    ")"
    ]
   },
   {
@@ -204,9 +204,7 @@
     "    energy_bins, gamma_rate, background_rate, config, assumed_spectrum\n",
     ")\n",
     "\n",
-    "combined_significance = source_detection(\n",
-    "    sigmas, u.Quantity(config[\"observation_time\"])\n",
-    ")"
+    "combined_significance = source_detection(sigmas, u.Quantity(config[\"observation_time\"]))"
    ]
   },
   {
@@ -215,33 +213,42 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "annotation_options = {\"rotation\": 45,\n",
-    "        \"xytext\": (10, 10),\n",
-    "        \"size\": 15}\n",
+    "annotation_options = {\"rotation\": 45, \"xytext\": (10, 10), \"size\": 15}\n",
     "\n",
     "with quantity_support():\n",
     "    plot_sed(\n",
-    "            config,\n",
-    "            sigmas,\n",
-    "            combined_significance,\n",
-    "            source_name,\n",
-    "            assumed_spectrum,\n",
-    "            en,\n",
-    "            sed,\n",
-    "            dsed,\n",
-    "            detected,\n",
-    "            savefig=False,\n",
-    "            annotation_options=annotation_options,\n",
-    "            )\n",
-    "    plt.ylim(1.e-12, 2.e-10)"
+    "        config,\n",
+    "        sigmas,\n",
+    "        combined_significance,\n",
+    "        source_name,\n",
+    "        assumed_spectrum,\n",
+    "        en,\n",
+    "        sed,\n",
+    "        dsed,\n",
+    "        detected,\n",
+    "        savefig=False,\n",
+    "        annotation_options=annotation_options,\n",
+    "    )\n",
+    "    plt.ylim(1.0e-12, 2.0e-10)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "instrument_fov = u.Quantity(config[\"fov\"])\n",
+    "wobble_offsets, wobble_angles = load_wobbles(config[\"wobbles\"])\n",
+    "plot_skymap_with_wobbles(\n",
+    "    target_source,\n",
+    "    observer,\n",
+    "    instrument_fov,\n",
+    "    wobble_angles,\n",
+    "    wobble_offsets,\n",
+    "    config,\n",
+    ")"
+   ]
   }
  ],
  "metadata": {

environment.yml

@@ -30,3 +30,6 @@ dependencies:
   # Development
   - pre-commit
   - ruff
+  # PyPI-only or dev dependencies
+  - pip:
+      - starplot

pyproject.toml

@@ -28,6 +28,7 @@ dependencies = [
     "gammapy",
     "seaborn",
     "scipy<1.12", # see https://github.com/gammapy/gammapy/pull/4997
+    "starplot"
 ]
 dynamic = ["version"]
 

src/iact_estimator/io.py

@@ -4,12 +4,7 @@
 from pathlib import Path
 
 from numpy import loadtxt
-from yaml import load
-
-try:
-    from yaml import CLoader as Loader
-except ImportError:
-    from yaml import Loader
+from yaml import safe_load
 
 __all__ = ["read_yaml", "load_ebl"]
 
@@ -33,7 +28,7 @@ def read_yaml(input_file_path):
     """
     try:
         with open(input_file_path, "r") as input_file:
-            data = load(input_file, Loader=Loader)
+            data = safe_load(input_file)
     except FileNotFoundError:
         logger.exception("Configuration file not found at %s", input_file_path)
     return data

src/iact_estimator/plots/observability.py

@@ -0,0 +1,149 @@
+"""Plotting functions related to source observability from a location."""
+
+import logging
+from pathlib import Path
+
+import astropy.units as u
+from astroplan import FixedTarget
+from astroplan.plots import plot_sky_24hr, plot_altitude
+import matplotlib.pyplot as plt
+
+from ..core import get_horizon_stereo_profile
+from iact_estimator import HORIZON_PROFILE_M1, HORIZON_PROFILE_M2
+
+logger = logging.getLogger(__name__)
+
+__all__ = ["plot_transit", "plot_altitude_airmass"]
+
+
+def plot_transit(
+    config,
+    source_name,
+    target_source,
+    observer,
+    time,
+    merge_profiles=True,
+    plot_crab=True,
+    style_kwargs=None,
+    savefig=True,
+    output_path=None,
+):
+    """
+    Plot the Spectral Energy distribution with significances.
+
+    Parameters
+    ----------
+    config : dict
+        Loaded configuration
+    output_path : str or `~pathlib.Path`
+    source_name : str
+    target_source : `astroplan.Target`
+    observer : `astroplan.Observer`
+    time : `~astropy.time.Time`
+        Datetime of the planned observation
+    merge_profiles : bool, default=True
+        If True plot the combined horizon profile
+        from both telescopes.
+    crab : bool, default=True
+        If True plot the Crab together with
+        the target source for comparison.
+    style_kwargs : dict, default=None
+        Dictionary of keywords passed into
+        `~matplotlib.pyplot.scatter`
+        to set plotting styles.
+    """
+    fig = plt.figure(figsize=config["plotting_options"]["figure_size"])
+    ax = fig.add_subplot(projection="polar")
+    ax.tick_params(pad=10)
+
+    plot_sky_24hr(
+        target_source,
+        observer,
+        time,
+        delta=1 * u.h,
+        ax=ax,
+        style_kwargs=style_kwargs,
+        north_to_east_ccw=True,
+        grid=True,
+        az_label_offset=0 * u.deg,
+        center_time_style_kwargs=None,
+    )
+    if plot_crab:
+        crab_nebula = FixedTarget.from_name("Crab Nebula")
+        plot_sky_24hr(
+            crab_nebula,
+            observer,
+            time,
+            delta=1 * u.h,
+            ax=ax,
+            style_kwargs={"label": "Crab Nebula"},
+            north_to_east_ccw=True,
+            grid=True,
+            az_label_offset=0 * u.deg,
+            center_time_style_kwargs=None,
+        )
+
+    if merge_profiles:
+        az, zd = get_horizon_stereo_profile(HORIZON_PROFILE_M1, HORIZON_PROFILE_M2)
+        ax.plot(
+            az.to_value("rad"),
+            zd.to_value("deg"),
+            linestyle="-",
+            label="Horizon profile",
+            alpha=0.5,
+            color="#4daf4a",
+        )
+    else:
+        ax.plot(
+            HORIZON_PROFILE_M2["azimuth"].to("rad").value,
+            HORIZON_PROFILE_M2["zenith"].value,
+            linestyle="-",
+            label="Horizon profile M2",
+            alpha=0.5,
+            color="#4daf4a",
+        )
+        ax.plot(
+            HORIZON_PROFILE_M1["azimuth"].to("rad").value,
+            HORIZON_PROFILE_M1["zenith"].value,
+            linestyle="-",
+            label="Horizon profile M1",
+            alpha=0.5,
+            color="#f781bf",
+        )
+
+    ax.legend(loc="center", bbox_to_anchor=(0.5, 0.8))
+    if savefig:
+        output_path = output_path if output_path is not None else Path.cwd()
+        fig.savefig(
+            output_path
+            / f"{source_name}_skyplot.{config['plotting_options']['file_format']}",
+            bbox_inches=config["plotting_options"]["bbox_inches"],
+        )
+        logger.debug("Plot has been successfully saved at %s", output_path)
+    return ax
+
+
+def plot_altitude_airmass(
+    config,
+    source_name,
+    target_source,
+    observer,
+    time,
+    ax=None,
+    savefig=True,
+    output_path=None,
+    **kwargs,
+):
+    fig, ax = plt.subplots(figsize=config["plotting_options"]["figure_size"])
+
+    plot_altitude(target_source, observer, time, ax=ax, **kwargs)
+    ax.grid(False)
+
+    if savefig:
+        output_path = output_path if output_path is not None else Path.cwd()
+        fig.savefig(
+            output_path
+            / f"{source_name}_altitude_airmass.{config['plotting_options']['file_format']}",
+            bbox_inches=config["plotting_options"]["bbox_inches"],
+        )
+    return ax