diff --git a/docs/source/changes/25.added.rst b/docs/source/changes/25.added.rst new file mode 100644 index 0000000..7dd22f8 --- /dev/null +++ b/docs/source/changes/25.added.rst @@ -0,0 +1 @@ +Add support for plotting and saving images from the `SkyView `_ service. diff --git a/docs/source/userguide.rst b/docs/source/userguide.rst index 5a189f9..e346417 100644 --- a/docs/source/userguide.rst +++ b/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 ====== diff --git a/environment.yml b/environment.yml index 15917bd..ffc4182 100644 --- a/environment.yml +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index edd2e80..a59dfd6 100644 --- a/pyproject.toml +++ b/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 diff --git a/src/iact_estimator/io.py b/src/iact_estimator/io.py index f9a790d..4f1fa40 100644 --- a/src/iact_estimator/io.py +++ b/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) diff --git a/src/iact_estimator/plots.py b/src/iact_estimator/plots.py index 8af48a9..7a7e60d 100644 --- a/src/iact_estimator/plots.py +++ b/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 diff --git a/src/iact_estimator/resources/config.yml b/src/iact_estimator/resources/config.yml index 45190c4..695dd0e 100644 --- a/src/iact_estimator/resources/config.yml +++ b/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 diff --git a/src/iact_estimator/scripts/main.py b/src/iact_estimator/scripts/main.py index e9cb4a5..cc4be2b 100644 --- a/src/iact_estimator/scripts/main.py +++ b/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)