Add checks for observability and grid plot

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

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

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()