"""Plotting functions related to source observability from a location."""
import logging
import warnings
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from datetime import timedelta
from pathlib import Path
from astroplan.plots import plot_sky_24hr, plot_altitude
from astroplan.utils import time_grid_from_range
from astropy.coordinates.errors import NonRotationTransformationWarning
from astropy.time import Time
from ..observability import get_total_available_time
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"]
[docs]
def plot_transit(
config,
source_name,
target_source,
observer,
time,
merge_profiles=True,
style_kwargs=None,
savefig=True,
output_path=None,
):
"""
Plot the 24h skyplot of a target source as seen by an observer.
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 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
[docs]
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
def plot_observability_constraints_grid(
source_name,
config,
observer,
target,
start_time,
end_time,
time_resolution,
constraints,
ax=None,
savefig=True,
output_path=None,
):
"""
Plot a grid representing the observability of a target based on a set of constraints
over a specified time range.
Parameters
----------
source_name : str
The name of the source being observed.
config : dict
Configuration dictionary containing plotting options, such as file format.
observer : `~astroplan.Observer`
The observer object that provides the position and time of observation.
target : `~astroplanb.FixedTarget`
The target object representing the source of interest.
start_time : datetime
The starting time of the observation period.
end_time : datetime
The ending time of the observation period.
time_resolution : `~astropy.units.Quantity`
The time resolution (step size) for the grid.
constraints : `~astroplan.constraints.Constraint` or list(~astroplan.constraints.Constraint)`
A list of constraint functions that evaluate the observability of the target
at given times for the observer and target.
ax : `~matplotlib.axes.Axes`, optional
The axis on which to plot the grid. If None, a new figure and axis are created.
savefig : bool, optional
Whether to save the generated plot as a file. Default is True.
output_path : str or Path, optional
The directory where the plot should be saved. If None, saves to the current
working directory.
Returns
-------
ax : `~matplotlib.axes.Axes`
The axis with the generated observability grid plot.
Notes
-----
The observability grid is created by evaluating each constraint function
at every time step within the specified time range. The grid is plotted as
an image, with time on the x-axis and constraints on the y-axis.
"""
fig, ax = plt.subplots(layout="constrained")
ax = plt.cga() if ax is None else ax
# Create grid of times from ``start_time`` to ``end_time``
# with resolution ``time_resolution``
time_grid = time_grid_from_range(
[start_time, end_time], time_resolution=time_resolution
)
observability_grid = np.zeros((len(constraints), len(time_grid)))
constraint_labels = {}
for i, constraint in enumerate(constraints):
# Evaluate each constraint
observability_grid[i, :] = constraint(observer, target, times=time_grid)
constraint_labels[i] = f"{constraint.__class__.__name__.split('Constraint')[0]}"
# Create plot showing observability of the target
numcols = len(time_grid)
numrows = len(constraint_labels)
extent = [-0.5, numcols - 0.5, numrows - 0.5, -0.5]
ax.imshow(observability_grid, extent=extent)
ax.set_yticks(range(0, 3))
ax.set_yticks(list(constraint_labels.keys()), constraint_labels.values())
ax.set_xticks(range(len(time_grid)))
ax.set_xticklabels([t.datetime.strftime("%H:%M") for t in time_grid], fontsize=7)
ax.set_xticks(np.arange(extent[0], extent[1]), minor=True)
ax.set_yticks(np.arange(extent[2], extent[3]), minor=True)
ax.grid(which="minor", color="w", linestyle="-", linewidth=2)
ax.tick_params(axis="x", which="minor", bottom="off")
plt.setp(ax.get_xticklabels(), rotation=30, ha="right")
ax.tick_params(axis="y", which="minor", left="off")
if savefig:
output_path = output_path if output_path is not None else Path.cwd()
fig.savefig(
output_path
/ f"observability_grid_{source_name}.{config['plotting_options']['file_format']}",
)
return ax
def create_observability_heatmap(
target_source,
observer,
constraints,
start_date,
end_date,
time_resolution=1 * u.hour,
cmap="YlGnBu",
sns_plotting_context="paper",
sns_axes_style="whitegrid",
savefig=True,
output_path=None,
save_format="png",
):
"""Plot an annotated heatmap showing the amount of available hours per day.
Parameters
==========
target_source: `~astroplan.FixedTarget`
observer: `~astroplan.Observer`
constraints: `~astroplan.Constraint` or `list(~astroplan.Constraint)`
start_date: `~astropy.time.Time`
end_date: `~astropy.time.Time`
time_resolution: `u.Quantity`, default=1h
cmap: str, default="YlGnBu"
sns_plotting_context: str, default="paper"
sns_axes_style: str, default="darkgrid"
savefig: bool, default=True
output_path: `pathlib.Path`, default=None
If unspecified the figure is saved in the current working directory.
save_format: str, default="png"
"""
date_range = pd.date_range(
start=start_date.datetime, end=end_date.datetime, freq="D"
)
data = []
for date in date_range:
# Define time range for this particular day
time_range = Time([date, date + timedelta(days=1)])
# Get available observation hours for the day
# see also https://github.com/astropy/astroplan/issues/598
warnings.filterwarnings("ignore", category=NonRotationTransformationWarning)
# see also https://github.com/astropy/astroplan/issues/132#issuecomment-225471962
warnings.filterwarnings(
"ignore", module="astropy.coordinates.builtin_frames.utils"
)
available_hours = get_total_available_time(
target_source, observer, constraints, time_range, time_resolution
).value
# Append data with month, day, and available hours
data.append(
{
"year": date.year,
"month": date.month,
"day": date.day,
"available_hours": available_hours,
}
)
# Create DataFrame
df = pd.DataFrame(data)
if df["year"].nunique() > 1:
df["month_label"] = (
df["year"].astype(str) + "-" + df["month"].astype(str).str.zfill(2)
)
else:
df["month_label"] = df["month"]
heatmap_data = df.pivot(
index="month_label", columns="day", values="available_hours"
)
with (
sns.axes_style(sns_axes_style),
sns.plotting_context(sns_plotting_context),
sns.color_palette("deep"),
):
# Dynamically set the figure size based on data dimensions
fig_width = 1 + heatmap_data.shape[1] * 0.3 # 0.3 inch per day
fig_height = 1 + heatmap_data.shape[0] * 0.5 # 0.5 inch per month
fig, ax = plt.subplots(figsize=(fig_width, fig_height), layout="constrained")
fmt = ".0f" if time_resolution == 1 * u.h else ".1f"
sns.heatmap(
heatmap_data,
ax=ax,
annot=True,
fmt=fmt,
cmap=cmap,
annot_kws={"size": 12, "fontweight": 10},
cbar_kws={"label": "Available Observation Hours"},
)
ax.set_xlabel("Day of the Month")
ax.set_ylabel("Month")
if savefig:
output_path = output_path if output_path is not None else Path.cwd()
fig.savefig(
output_path / f"observability_heatmap_{target_source.name}.{save_format}",
)