diff --git a/docs/source/changes/22.added.rst b/docs/source/changes/22.added.rst new file mode 100644 index 0000000..943c30c --- /dev/null +++ b/docs/source/changes/22.added.rst @@ -0,0 +1 @@ +Added preliminary checks for observability and new plot to show an observability grid in time based on the specified constraints. diff --git a/docs/source/example.ipynb b/docs/source/example.ipynb index 5b83a3a..2dbe897 100644 --- a/docs/source/example.ipynb +++ b/docs/source/example.ipynb @@ -29,13 +29,19 @@ "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.observability import define_constraints, get_days_in_this_year\n", + "from iact_estimator.plots import (\n", + " plot_spectrum,\n", + " plot_sed,\n", + " plot_transit,\n", + " plot_altitude_airmass,\n", + " create_observability_heatmap,\n", + ")" ] }, { @@ -46,8 +52,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 +103,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 +142,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 +179,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", + ")" ] }, { @@ -205,7 +210,7 @@ ")\n", "\n", "combined_significance = source_detection(\n", - " sigmas, u.Quantity(config[\"observation_time\"])\n", + " sigmas, u.Quantity(config[\"observation\"][\"time\"])\n", ")" ] }, @@ -215,25 +220,23 @@ "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)" ] }, { @@ -241,7 +244,36 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "constraints = define_constraints(config)\n", + "\n", + "start_datetime = (\n", + " Time(config[\"observation\"][\"start_datetime\"])\n", + " if config[\"observation\"][\"start_datetime\"] is not None\n", + " else Time(datetime.now(tz=observer.timezone))\n", + ")\n", + "year_days = get_days_in_this_year()\n", + "end_datetime = (\n", + " Time(config[\"observation\"][\"end_datetime\"])\n", + " if config[\"observation\"][\"end_datetime\"] is not None\n", + " else start_datetime + year_days\n", + ")\n", + "\n", + "create_observability_heatmap(\n", + " target_source,\n", + " observer,\n", + " constraints,\n", + " start_datetime,\n", + " end_datetime,\n", + " time_resolution=1 * u.hour,\n", + " cmap=\"YlGnBu\",\n", + " sns_plotting_context=\"paper\",\n", + " sns_axes_style=\"whitegrid\",\n", + " savefig=True,\n", + " output_path=None,\n", + " save_format=\"png\",\n", + ")" + ] } ], "metadata": { diff --git a/src/iact_estimator/core.py b/src/iact_estimator/core.py index c42b5ff..101ff43 100644 --- a/src/iact_estimator/core.py +++ b/src/iact_estimator/core.py @@ -429,7 +429,7 @@ def calculate(energy_bins, gamma_rate, background_rate, config, assumed_spectrum n_off_regions = config["n_off_regions"] redshift = config["redshift"] - observation_time_min = u.Quantity(config["observation_time"]).to("min") + observation_time_min = u.Quantity(config["observation"]["time"]).to("min") pulsar_mode_config = config["pulsar_mode"] pulsar_mode = pulsar_mode_config["enable"] pulsar_on_range = pulsar_mode_config["pulsar_on_range"] diff --git a/src/iact_estimator/observability.py b/src/iact_estimator/observability.py new file mode 100644 index 0000000..d7362a1 --- /dev/null +++ b/src/iact_estimator/observability.py @@ -0,0 +1,119 @@ +"""Module that stores all functions related to observability estimation from the given location.""" + +from datetime import datetime + +from astropy.time import Time + +from astroplan import ( + AtNightConstraint, + MoonSeparationConstraint, + AltitudeConstraint, + MoonIlluminationConstraint, +) +import astropy.units as u + +from astroplan import ( + is_observable, + months_observable, + is_event_observable, + time_grid_from_range, +) + + +__all__ = ["define_constraints"] + + +def get_days_in_this_year(): + """Get the number of days in the current year. + + Accounts for leap years.""" + this_year = datetime.today().year + start_of_year = Time(f"{this_year}-01-01") + end_of_year = Time(f"{this_year}-12-31") + + # Add 1 to include the last day + num_days = int((end_of_year - start_of_year).to_value("day") + 1) + return num_days + + +def define_constraints(config): + constraints_config = config["observation"]["constraints"] + + moon_separation = constraints_config["moon_separation"] + moon_illumination = constraints_config["moon_illumination_fraction"] + zenith = constraints_config["zenith"] + + constraints = [ + AtNightConstraint(u.Quantity(constraints_config["max_solar_altitude"])), + MoonSeparationConstraint( + min=u.Quantity(moon_separation["min"]), + max=u.Quantity(moon_separation["max"]) if moon_separation["max"] else None, + ephemeris=moon_separation["ephemeris"], + ), + MoonIlluminationConstraint( + min=u.Quantity(moon_illumination["min"]) + if moon_illumination["min"] + else None, + max=u.Quantity(moon_illumination["max"]), + ephemeris=moon_separation["ephemeris"], + ), + AltitudeConstraint( + min=90 * u.deg - u.Quantity(zenith["max"]), + max=90 * u.deg - u.Quantity(zenith["min"]), + ), + ] + + return constraints + + +def check_observability( + constraints, observer, target_source, time_range, time_grid_resolution +): + # Are targets *ever* observable in the time range? + ever_observable = is_observable( + constraints, + observer, + target_source, + time_range=time_range, + time_grid_resolution=time_grid_resolution, + ) + + # During what months are the targets ever observable? + best_months = months_observable( + constraints, + observer, + target_source, + time_range, + time_grid_resolution=time_grid_resolution, + ) + + return ever_observable, best_months + + +def get_total_available_time( + target_source, observer, constraints, time_range, time_resolution +): + """Get available observation time. + + Args: + target_source (astroplan.FixedTarget) + observer (astroplan.Observer) + constraints (astroplan.Constraint or list(astroplan.Constraint)): + Set of constraining limits for the observation. + time_range (list(astropy.time.Time)): [start, stop] time stamps + time_resolution (u.Quantity): time step over which to evaluate the constraints. + + Returns: + total_available_time (u.Quantity): in hours. + """ + time_grid = time_grid_from_range(time_range, time_resolution=time_resolution) + + is_observable_at_time = is_event_observable( + constraints, observer, target_source, times=time_grid + ) + + total_observable_time = is_observable_at_time.sum() + + total_available_time = total_observable_time * time_resolution.to("h") + + return total_available_time diff --git a/src/iact_estimator/plots.py b/src/iact_estimator/plots.py index 8af48a9..f5ffc7e 100644 --- a/src/iact_estimator/plots.py +++ b/src/iact_estimator/plots.py @@ -7,12 +7,18 @@ from astropy.visualization import quantity_support from astroplan import FixedTarget from astroplan.plots import plot_sky_24hr, plot_altitude +from astropy.time import Time +from datetime import timedelta +from astroplan.utils import time_grid_from_range import matplotlib.pyplot as plt import numpy as np +import pandas as pd +import seaborn as sns from .core import observed_flux, get_horizon_stereo_profile from .spectral import crab_nebula_spectrum from iact_estimator import HORIZON_PROFILE_M1, HORIZON_PROFILE_M2 +from .observability import get_total_available_time __all__ = [ "plot_spectrum", @@ -99,7 +105,7 @@ def plot_sed( "verticalalignment": "bottom", }, ): - """ + r""" Plot the Spectral Energy distribution with significances. Parameters @@ -156,7 +162,7 @@ def plot_sed( np.log10(max_energy.to_value(energy_unit)), 50, ) * u.Unit(energy_unit) - labeltext = rf"Expected SED ($T_{{obs}}$ = {config['observation_time']})" + labeltext = rf"Expected SED ($T_{{obs}}$ = {config['observation']['time']})" plt.plot( energy, energy * energy * crab_nebula_spectrum()(energy), @@ -188,7 +194,7 @@ def plot_sed( for i in range(len(sigmas)): col = "0" if detected[i] else "0.75" ax.annotate( - f"{sigmas[i]:.1f}$\sigma$", + rf"{sigmas[i]:.1f}$\sigma$", (en[i], sed[i]), color=col, xycoords="data", @@ -375,3 +381,164 @@ def plot_rates(performance_data, title=None, ax=None): ax.set_xscale("log") ax.set_yscale("log") 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, +): + 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 + 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}", + ) diff --git a/src/iact_estimator/resources/config.yml b/src/iact_estimator/resources/config.yml index 45190c4..a7c969d 100644 --- a/src/iact_estimator/resources/config.yml +++ b/src/iact_estimator/resources/config.yml @@ -8,8 +8,26 @@ assumed_model: beta: 0.21 from_log10: True # relevant for e.g. LogParabolaSpectralModel -observation_time: 50 h -observation_datetime: "2024-06-15 18:00" +observation: + constraints: + moon_illumination_fraction: + max: 0.2 + min: + ephemeris: + moon_separation: + # https://astroplan.readthedocs.io/en/latest/api/astroplan.MoonSeparationConstraint.html#astroplan.MoonSeparationConstraint + min: "45 deg" + max: + ephemeris: + zenith: + min: "0 deg" + max: "60 deg" + max_solar_altitude: "-18 deg" + time: "50 h" + start_datetime: # default to today + end_datetime: # default to 1 year from today + time_resolution: "1 h" + extension: 0.0 deg redshift: -1 # redshift of the source (for the EBL absorption), if -1 then no absorption sum_trigger: False diff --git a/src/iact_estimator/scripts/main.py b/src/iact_estimator/scripts/main.py index e9cb4a5..8217ea1 100644 --- a/src/iact_estimator/scripts/main.py +++ b/src/iact_estimator/scripts/main.py @@ -3,6 +3,7 @@ import argparse import logging from pathlib import Path +from datetime import datetime import shutil from astroplan import FixedTarget, Observer @@ -21,7 +22,19 @@ 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_observability_constraints_grid, + create_observability_heatmap, +) +from ..observability import ( + define_constraints, + check_observability, + get_days_in_this_year, +) from .. import RESOURCES_PATH parser = argparse.ArgumentParser() @@ -127,6 +140,109 @@ def main(): seaborn_options = config["seaborn_options"] sns.set_theme(**seaborn_options) + # Basic observability checks + target_source = FixedTarget.from_name(source_name) + observer = Observer.at_site("Roque de los Muchachos") + + crab = FixedTarget.from_name("Crab") + + logger.debug("Defining observation constraints") + constraints = define_constraints(config) + + start_datetime = ( + Time(config["observation"]["start_datetime"]) + if config["observation"]["start_datetime"] is not None + else Time(datetime.now(tz=observer.timezone)) + ) + year_days = get_days_in_this_year() + end_datetime = ( + Time(config["observation"]["end_datetime"]) + if config["observation"]["end_datetime"] is not None + else start_datetime + year_days + ) + logger.info("Observation starts at %s", start_datetime) + logger.info("Observation ends at %s", end_datetime) + + time_range = [start_datetime, end_datetime] + time_grid_resolution = ( + u.Quantity(config["observation"]["time_resolution"]) + if config["observation"]["time_resolution"] + else 1 * u.h + ) + + logger.debug("Checking observability") + ever_observable, best_months = check_observability( + constraints, observer, [target_source], time_range, time_grid_resolution + ) + + logger.debug("Producing observability constraints grid") + obs_grid_time_res = ( + 1 * u.h + if (end_datetime - start_datetime).to("day").value <= 1 + else 1 * u.day + ) + _ = plot_observability_constraints_grid( + source_name, + config, + observer, + target_source, + start_datetime, + end_datetime, + obs_grid_time_res, + constraints, + ax=None, + savefig=True, + output_path=output_path, + ) + + if not ever_observable: + logger.info("The source is never observable from this location!") + parser.exit(0) + + logger.info(f"The best months to observe the target source are {best_months}.") + + logger.debug("Producing observability heatmap") + create_observability_heatmap( + target_source, + observer, + constraints, + start_datetime, + end_datetime, + time_resolution=1 * u.hour, + cmap="YlGnBu", + sns_plotting_context="paper", + sns_axes_style="whitegrid", + savefig=True, + output_path=None, + save_format="png", + ) + + with quantity_support(): + plot_transit( + config, + source_name, + target_source, + observer, + start_datetime, + merge_profiles=config["plotting_options"]["merge_horizon_profiles"], + plot_crab=True if (crab.coord == target_source.coord) else False, + style_kwargs=None, + savefig=True, + output_path=output_path, + ) + + plot_altitude_airmass( + config, + source_name, + target_source, + observer, + start_datetime, + brightness_shading=True, + airmass_yaxis=True, + savefig=True, + output_path=output_path, + ) + logger.info("Initializing assumed model") assumed_spectrum = initialize_model(config) @@ -155,7 +271,7 @@ def main(): ) combined_significance = source_detection( - sigmas, u.Quantity(config["observation_time"]) + sigmas, u.Quantity(config["observation"]["time"]) ) with quantity_support(): @@ -175,38 +291,6 @@ def main(): logger.info("All expected operations have been perfomed succesfully.") - target_source = FixedTarget.from_name(source_name) - observer = Observer.at_site("Roque de los Muchachos") - time = Time(config["observation_datetime"]) - - crab = FixedTarget.from_name("Crab") - - with quantity_support(): - plot_transit( - config, - source_name, - target_source, - observer, - time, - merge_profiles=config["plotting_options"]["merge_horizon_profiles"], - plot_crab=True if (crab.coord == target_source.coord) else False, - style_kwargs=None, - savefig=True, - output_path=output_path, - ) - - plot_altitude_airmass( - config, - source_name, - target_source, - observer, - time, - brightness_shading=True, - airmass_yaxis=True, - savefig=True, - output_path=output_path, - ) - if config["plotting_options"]["show"]: plt.show()