From 14778b7ef4b97e7f71bf06da0f61492e09f91c96 Mon Sep 17 00:00:00 2001 From: Michele Peresano Date: Thu, 19 Sep 2024 17:18:20 +0200 Subject: [PATCH 1/3] Add new observability module --- src/iact_estimator/observability.py | 41 +++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 src/iact_estimator/observability.py diff --git a/src/iact_estimator/observability.py b/src/iact_estimator/observability.py new file mode 100644 index 0000000..99cb80a --- /dev/null +++ b/src/iact_estimator/observability.py @@ -0,0 +1,41 @@ +"""Module that stores all functions related to observability estimation from the given location.""" + +from astroplan import AtNightConstraint, MoonSeparationConstraint, AltitudeConstraint +import astropy.units as u + +from astroplan import is_observable, is_always_observable, months_observable + + +__all__ = ["define_constraints"] + + +def define_constraints(config): + constraints_config = config["observation"]["constraints"] + + constraints = [ + getattr(AtNightConstraint, f"twilight_{constraints_config["twilight"]}"), + MoonSeparationConstraint(**constraints_config["moon_separation"]), + AltitudeConstraint( + min=90 * u.deg - constraints_config["zenith"]["max"], + max=90 * u.deg - constraints_config["zenith"]["min"], + ), + ] + + return constraints + + +def check_observability(constraints, observer, target_source, time_range): + # Are targets *ever* observable in the time range? + ever_observable = is_observable( + constraints, observer, target_source, time_range=time_range + ) + + # Are targets *always* observable in the time range? + always_observable = is_always_observable( + constraints, observer, target_source, time_range=time_range + ) + + # During what months are the targets ever observable? + best_months = months_observable(constraints, observer, target_source, time_range) + + return ever_observable, always_observable, best_months From a5f8252346a0045e721a61ed9eef0f2f7d861aac Mon Sep 17 00:00:00 2001 From: Michele Peresano Date: Thu, 19 Sep 2024 17:18:59 +0200 Subject: [PATCH 2/3] Add plot_observability_constraints_grid --- src/iact_estimator/plots.py | 43 +++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/src/iact_estimator/plots.py b/src/iact_estimator/plots.py index 8af48a9..37ff846 100644 --- a/src/iact_estimator/plots.py +++ b/src/iact_estimator/plots.py @@ -7,6 +7,7 @@ from astropy.visualization import quantity_support from astroplan import FixedTarget from astroplan.plots import plot_sky_24hr, plot_altitude +from astroplan.utils import time_grid_from_range import matplotlib.pyplot as plt import numpy as np @@ -375,3 +376,45 @@ def plot_rates(performance_data, title=None, ax=None): ax.set_xscale("log") ax.set_yscale("log") return ax + + +def plot_observability_constraints_grid( + observer, target, start_time, end_time, time_resolution, constraints, ax=None +): + # 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))) + + for i, constraint in enumerate(constraints): + # Evaluate each constraint + observability_grid[i, :] = constraint(observer, target, times=time_grid) + + # Create plot showing observability of the target + + extent = [-0.5, -0.5 + len(time_grid), -0.5, 2.5] + + ax = plt.cga() if ax is None else ax + + ax.imshow(observability_grid, extent=extent) + + ax.set_yticks(range(0, 3)) + ax.set_yticklabels([c.__class__.__name__ for c in constraints]) + + ax.set_xticks(range(len(time_grid))) + ax.set_xticklabels([t.datetime.strftime("%H:%M") for t in time_grid]) + + 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") + ax.set_xlabel("Time on {0} UTC".format(time_grid[0].datetime.date())) + + return ax From 1f6da457c5714d1161780d2b21b1b293585d1dda Mon Sep 17 00:00:00 2001 From: Michele Peresano Date: Thu, 19 Sep 2024 17:19:11 +0200 Subject: [PATCH 3/3] Update main with observability check --- src/iact_estimator/scripts/main.py | 80 ++++++++++++++++++------------ 1 file changed, 48 insertions(+), 32 deletions(-) diff --git a/src/iact_estimator/scripts/main.py b/src/iact_estimator/scripts/main.py index e9cb4a5..30b5035 100644 --- a/src/iact_estimator/scripts/main.py +++ b/src/iact_estimator/scripts/main.py @@ -22,6 +22,7 @@ calculate, ) from ..plots import plot_spectrum, plot_sed, plot_transit, plot_altitude_airmass +from ..observability import define_constraints, check_observability from .. import RESOURCES_PATH parser = argparse.ArgumentParser() @@ -127,6 +128,53 @@ 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") + time = Time(config["observation_datetime"]) + + crab = FixedTarget.from_name("Crab") + + constraints = define_constraints(config) + + time_range = config["observation"]["datetime"] + 1 * u.year + + ever_observable, always_observable, best_months = check_observability( + constraints, observer, target_source, time_range + ) + + 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}.") + + 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, + ) + logger.info("Initializing assumed model") assumed_spectrum = initialize_model(config) @@ -175,38 +223,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()