Add support for SkyView service

docs/source/changes/25.added.rst

@@ -0,0 +1 @@
+Add support for plotting and saving images from the `SkyView <https://skyview.gsfc.nasa.gov/current/cgi/titlepage.pl>`_ service.

docs/source/userguide.rst

@@ -11,7 +11,7 @@ You can use the :ref:`iact-estimator-cfg` command-line tool
 
 .. code-block::
 
-    iact-estimator-cfg --output-path /where/to/save/config/file
+    iact-estimator config --to /where/to/save/config/file
 
 to get the following example configuration file,
 
@@ -30,7 +30,7 @@ for a source named "my_source",
 
 .. code-block::
 
-    iact-estimator --config config.yml --source-name my_source
+    iact-estimator run --config config.yml --source-name my_source
 
 Output
 ======

environment.yml

@@ -9,6 +9,7 @@ dependencies:
   - pip
   - astroplan
   - astropy<6 # see https://github.com/gammapy/gammapy/issues/4972
+  - astroquery
   - gammapy
   - seaborn
   - scipy=1.11* # see https://github.com/gammapy/gammapy/pull/4997

pyproject.toml

@@ -24,6 +24,7 @@ classifiers = [
 dependencies = [
     "astroplan",
     "astropy<6",  # see https://github.com/gammapy/gammapy/issues/4972
+    "astroquery",
     "gammapy",
     "seaborn",
     "scipy<1.12", # see https://github.com/gammapy/gammapy/pull/4997

src/iact_estimator/io.py

@@ -75,3 +75,23 @@ def load_ebl(ebl_file_path):
     taus = ebl_body[:, 1:]
     if len(taus > 0):
         return zz, energies, taus
+
+
+def save_fits_hdu(hdu, output_path, **kwargs):
+    """Save and HDU to a new FITS file.
+
+    Parameters
+    ----------
+    hdu : HDU-like class from `astropy.io.fits`
+        Primary, Image or Table-HDU.
+    output_path : str or `~pathlib.Path`
+        Complete path to the saved file.
+    **kwargs : dict
+        Options for the `~astropy.io.fits.PrimaryHDU.writeto()` method.
+    """
+
+    if not Path(output_path).parent.is_dir():
+        raise ValueError("Output directory for %s does not exist.", output_path)
+
+    logger.info("Saving HDU to %s", output_path)
+    hdu.writeto(output_path, **kwargs)

src/iact_estimator/plots.py

@@ -6,7 +6,7 @@
 import astropy.units as u
 from astropy.visualization import quantity_support
 from astroplan import FixedTarget
-from astroplan.plots import plot_sky_24hr, plot_altitude
+from astroplan.plots import plot_sky_24hr, plot_altitude, plot_finder_image
 import matplotlib.pyplot as plt
 import numpy as np
 
@@ -375,3 +375,84 @@ def plot_rates(performance_data, title=None, ax=None):
         ax.set_xscale("log")
         ax.set_yscale("log")
     return ax
+
+
+@u.quantity_input(fov_radius=u.arcmin)
+def plot_from_skyview_survey(
+    target_source,
+    survey_name="DSS",
+    fov_radius=10 * u.arcmin,
+    log=False,
+    ax=None,
+    grid=False,
+    reticle=False,
+    style_kwargs=None,
+    reticle_style_kwargs=None,
+):
+    """Plot a survey image from the SkyView service centered on ``target``.
+
+    Parameters
+    ----------
+    target : `~astroplan.FixedTarget`, `~astropy.coordinates.SkyCoord`
+        Coordinates of celestial object
+
+    survey : string
+        Name of survey to retrieve image from. For dictionary of
+        available surveys, use
+        ``from astroquery.skyview import SkyView; SkyView.list_surveys()``.
+        Defaults to ``'DSS'``, the Digital Sky Survey.
+
+    fov_radius : `~astropy.units.Quantity`
+        Radius of field of view of retrieved image. Defaults to 10 arcmin.
+
+    log : bool, optional
+        Take the natural logarithm of the FITS image if `True`.
+        False by default.
+
+    ax : `~matplotlib.axes.Axes` or None, optional.
+        The `~matplotlib.axes.Axes` object to be drawn on.
+        If None, uses the current `~matplotlib.axes.Axes`.
+
+    grid : bool, optional.
+        Grid is drawn if `True`. `False` by default.
+
+    reticle : bool, optional
+        Draw reticle on the center of the FOV if `True`. Default is `False`.
+
+    style_kwargs : dict or `None`, optional.
+        A dictionary of keywords passed into `~matplotlib.pyplot.imshow`
+        to set plotting styles.
+
+    reticle_style_kwargs : dict or `None`, optional
+        A dictionary of keywords passed into `~matplotlib.pyplot.axvline` and
+        `~matplotlib.pyplot.axhline` to set reticle style.
+
+    Returns
+    -------
+    ax : `~matplotlib.axes.Axes`
+        Matplotlib axes with survey image centered on ``target``
+
+    hdu : `~astropy.io.fits.PrimaryHDU`
+        FITS HDU of the retrieved image
+
+    Notes
+    -----
+    This is wrapper function around `astroplan.plots.plot_finder_image()`.
+
+    """
+
+    ax = plt.gca() if ax is None else ax
+
+    ax = plot_finder_image(
+        target_source,
+        survey=survey_name,
+        fov_radius=fov_radius,
+        log=log,
+        ax=ax,
+        grid=grid,
+        reticle=reticle,
+        style_kwargs=style_kwargs,
+        reticle_style_kwargs=reticle_style_kwargs,
+    )
+
+    return ax

src/iact_estimator/resources/config.yml

@@ -44,6 +44,24 @@ plotting_options:
   file_format: "pdf"
   merge_horizon_profiles: True
 
+skyview:
+  save_hdus: True
+  surveys:
+    - name: "DSS"
+      fov_radius: "10 arcmin"
+      log: False
+      grid: False
+      reticle: False
+      style_kwargs: {"cmap": "viridis"}
+      reticle_style_kwargs: {"color": "white", "lw":2}
+    - name: "GB6 (4850MHz)"
+      fov_radius: "10 arcmin"
+      log: False
+      grid: False
+      reticle: False
+      style_kwargs: {"cmap": "viridis"}
+      reticle_style_kwargs: {"color": "white", "lw":2}
+
 use_seaborn: True
 seaborn_options:
   context: talk

src/iact_estimator/scripts/main.py

@@ -12,7 +12,7 @@
 import matplotlib.pyplot as plt
 
 from .. import __version__
-from ..io import read_yaml
+from ..io import read_yaml, save_fits_hdu
 from ..core import (
     setup_logging,
     initialize_model,
@@ -21,7 +21,13 @@
     source_detection,
     calculate,
 )
-from ..plots import plot_spectrum, plot_sed, plot_transit, plot_altitude_airmass
+from ..plots import (
+    plot_spectrum,
+    plot_sed,
+    plot_transit,
+    plot_altitude_airmass,
+    plot_from_skyview_survey,
+)
 from .. import RESOURCES_PATH
 
 parser = argparse.ArgumentParser()
@@ -127,6 +133,36 @@ def main():
             seaborn_options = config["seaborn_options"]
             sns.set_theme(**seaborn_options)
 
+        target_source = FixedTarget.from_name(source_name)
+
+        for survey in config["skyview"]["surveys"]:
+            survey_name = survey["name"]
+            fig, ax = plt.subplots()
+            _, hdu = plot_from_skyview_survey(
+                target_source,
+                survey_name=survey_name,
+                fov_radius=u.Quantity(survey["fov_radius"]),
+                log=survey["log"],
+                ax=ax,
+                grid=survey["grid"],
+                reticle=survey["reticle"],
+                style_kwargs=survey["style_kwargs"],
+                reticle_style_kwargs=survey["reticle_style_kwargs"],
+            )
+            ax.set_title(f"{source_name} - {survey_name}")
+            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"],
+            )
+            if config["skyview"]["save_hdus"]:
+                save_fits_hdu(
+                    hdu,
+                    output_path
+                    / f"{source_name}_skyview_image_from_{survey_name.replace(' ','')}.fits",
+                )
+
         logger.info("Initializing assumed model")
         assumed_spectrum = initialize_model(config)