Source code for astroplan.constraints

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
Specify and constraints to determine which targets are observable for
an observer.
"""

from __future__ import (absolute_import, division, print_function,
                        unicode_literals)

# Standard library
from abc import ABCMeta, abstractmethod
import datetime

# Third-party
from astropy.time import Time
import astropy.units as u
from astropy.coordinates import get_sun, Angle, SkyCoord
from astropy.coordinates.angle_utilities import angular_separation
from astropy import table
import numpy as np

# Package
from .moon import get_moon, moon_illumination
from .utils import time_grid_from_range
from .target import FixedTarget

__all__ = ["AltitudeConstraint", "AirmassConstraint", "AtNightConstraint",
           "is_observable", "is_always_observable", "time_grid_from_range",
           "SunSeparationConstraint", "MoonSeparationConstraint",
           "MoonIlluminationConstraint", "LocalTimeConstraint", "Constraint",
           "observability_table"]


def _get_altaz(times, observer, targets,
               force_zero_pressure=False):
    """
    Calculate alt/az for ``target`` at times linearly spaced between
    the two times in ``time_range`` with grid spacing ``time_resolution``
    for ``observer``.

    Cache the result on the ``observer`` object.

    Parameters
    ----------
    times : `~astropy.time.Time`
        Array of times on which to test the constraint

    targets : {list, `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`}
        Target or list of targets

    observer : `~astroplan.Observer`
        The observer who has constraints ``constraints``

    time_resolution : `~astropy.units.Quantity` (optional)
        Set the time resolution in calculations of the altitude/azimuth

    Returns
    -------
    altaz_dict : dict
        Dictionary containing two key-value pairs. (1) 'times' contains the
        times for the alt/az computations, (2) 'altaz' contains the
        corresponding alt/az coordinates at those times.
    """
    if not hasattr(observer, '_altaz_cache'):
        observer._altaz_cache = {}

    # convert times, targets to tuple for hashing
    aakey = (tuple(times.jd), tuple(targets))

    if aakey not in observer._altaz_cache:
        try:
            if force_zero_pressure:
                observer_old_pressure = observer.pressure
                observer.pressure = 0

            altaz = observer.altaz(times, targets)
            observer._altaz_cache[aakey] = dict(times=times,
                                                altaz=altaz)
        finally:
            if force_zero_pressure:
                observer.pressure = observer_old_pressure

    return observer._altaz_cache[aakey]


@abstractmethod
[docs]class Constraint(object): """ Abstract class for objects defining observational constraints. """ __metaclass__ = ABCMeta
[docs] def __call__(self, observer, targets, times=None, time_range=None, time_grid_resolution=0.5*u.hour): if times is None and time_range is not None: times = time_grid_from_range(time_range, time_resolution=time_grid_resolution) elif not isinstance(times, Time): times = Time(times) if times.isscalar: times = Time([times]) if hasattr(targets, '__len__'): targets = [FixedTarget(coord=target) if isinstance(target, SkyCoord) else target for target in targets] else: if isinstance(targets, SkyCoord): targets = FixedTarget(coord=targets) cons = self.compute_constraint(times, observer, targets) return cons
[docs] def compute_constraint(self, times, observer, targets): # Should be implemented on each subclass of Constraint raise NotImplementedError
[docs]class AltitudeConstraint(Constraint): """ Constrain the altitude of the target. .. note:: This will misbehave if you try to constrain negative altitudes, as the `~astropy.coordinates.AltAz` frame tends to mishandle negative altitudes. """ def __init__(self, min=None, max=None): """ Parameters ---------- min : `~astropy.units.Quantity` or `None` Minimum altitude of the target. `None` indicates no limit. max : `~astropy.units.Quantity` or `None` Maximum altitude of the target. `None` indicates no limit. """ if min is None: self.min = 0*u.deg else: self.min = min if max is None: self.max = 90*u.deg else: self.max = max
[docs] def compute_constraint(self, times, observer, targets): cached_altaz = _get_altaz(times, observer, targets) altaz = cached_altaz['altaz'] lowermask = self.min < altaz.alt uppermask = altaz.alt < self.max return lowermask & uppermask
[docs]class AirmassConstraint(AltitudeConstraint): """ Constrain the airmass of a target. In the current implementation the airmass is approximated by the secant of the zenith angle. """ def __init__(self, max=None, min=None): """ .. note:: The ``max`` and ``min`` arguments appear in the order (max, min) in this initializer to support the common case for users who care about the upper limit on the airmass (``max``) and not the lower limit. Parameters ---------- max : float or `None` Maximum airmass of the target. `None` indicates no limit. min : float or `None` Minimum airmass of the target. `None` indicates no limit. Examples -------- To create a constraint that requires the airmass be "better than 2", i.e. at a higher altitude than airmass=2:: AirmassConstraint(2) """ self.min = min self.max = max
[docs] def compute_constraint(self, times, observer, targets): cached_altaz = _get_altaz(times, observer, targets) altaz = cached_altaz['altaz'] if self.min is None and self.max is not None: mask = altaz.secz < self.max elif self.max is None and self.min is not None: mask = self.min < altaz.secz elif self.min is not None and self.max is not None: mask = (self.min < altaz.secz) & (altaz.secz < self.max) else: raise ValueError("No max and/or min specified in " "AirmassConstraint.") return mask
[docs]class AtNightConstraint(Constraint): """ Constrain the Sun to be below ``horizon``. """ @u.quantity_input(horizon=u.deg) def __init__(self, max_solar_altitude=0*u.deg, force_pressure_zero=True): """ Parameters ---------- max_solar_altitude : `~astropy.units.Quantity` Define "night" as when the sun is below ``max_solar_altitude``. Default is zero degrees altitude. force_pressure_zero : bool (optional) Force the pressure to zero for solar altitude calculations. This avoids errors in the altitude of the Sun that can occur when the Sun is below the horizon and the corrections for atmospheric refraction return nonsense values. Default is `True`. """ self.max_solar_altitude = max_solar_altitude self.force_pressure_zero = force_pressure_zero @classmethod
[docs] def twilight_civil(cls, **kwargs): """ Consider nighttime as time between civil twilights (-6 degrees). """ return cls(max_solar_altitude=-6*u.deg, **kwargs)
@classmethod
[docs] def twilight_nautical(cls, **kwargs): """ Consider nighttime as time between nautical twilights (-12 degrees). """ return cls(max_solar_altitude=-12*u.deg, **kwargs)
@classmethod
[docs] def twilight_astronomical(cls, **kwargs): """ Consider nighttime as time between astronomical twilights (-18 degrees). """ return cls(max_solar_altitude=-18*u.deg, **kwargs)
def _get_solar_altitudes(self, times, observer, targets): if not hasattr(observer, '_altaz_cache'): observer._altaz_cache = {} aakey = (tuple(times.jd), 'sun') if aakey not in observer._altaz_cache: try: if self.force_pressure_zero: observer_old_pressure = observer.pressure observer.pressure = 0 # Broadcast the solar altitudes for the number of targets: altaz = observer.altaz(times, get_sun(times)) altitude = altaz.alt altitude.resize(1, len(altitude)) altitude = altitude + np.zeros((len(targets), 1)) observer._altaz_cache[aakey] = dict(times=times, altitude=altitude) finally: if self.force_pressure_zero: observer.pressure = observer_old_pressure return observer._altaz_cache[aakey]
[docs] def compute_constraint(self, times, observer, targets): sun_altaz = self._get_solar_altitudes(times, observer, targets) solar_altitude = sun_altaz['altitude'] mask = solar_altitude < self.max_solar_altitude return mask
[docs]class SunSeparationConstraint(Constraint): """ Constrain the distance between the Sun and some targets. """ def __init__(self, min=None, max=None): """ Parameters ---------- min : `~astropy.units.Quantity` or `None` (optional) Minimum acceptable separation between Sun and target. `None` indicates no limit. max : `~astropy.units.Quantity` or `None` (optional) Minimum acceptable separation between Sun and target. `None` indicates no limit. """ self.min = min self.max = max
[docs] def compute_constraint(self, times, observer, targets): sun = get_sun(times) targets = [target.coord if hasattr(target, 'coord') else target for target in targets] solar_separation = Angle([sun.separation(target) for target in targets]) if self.min is None and self.max is not None: mask = self.max > solar_separation elif self.max is None and self.min is not None: mask = self.min < solar_separation elif self.min is not None and self.max is not None: mask = ((self.min < solar_separation) & (solar_separation < self.max)) else: raise ValueError("No max and/or min specified in " "SunSeparationConstraint.") return mask
[docs]class MoonSeparationConstraint(Constraint): """ Constrain the distance between the Earth's moon and some targets. """ def __init__(self, min=None, max=None): """ Parameters ---------- min : `~astropy.units.Quantity` or `None` (optional) Minimum acceptable separation between moon and target. `None` indicates no limit. max : `~astropy.units.Quantity` or `None` (optional) Minimum acceptable separation between moon and target. `None` indicates no limit. """ self.min = min self.max = max
[docs] def compute_constraint(self, times, observer, targets): moon = get_moon(times, observer.location, observer.pressure) targets = [target.coord if hasattr(target, 'coord') else target for target in targets] moon_separation = Angle([angular_separation(moon.spherical.lon, moon.spherical.lat, target.spherical.lon, target.spherical.lat) for target in targets]) # The line below should have worked, but needs a workaround. # TODO: once bug has been fixed, replace workaround with simpler version. # Relevant PR: https://github.com/astropy/astropy/issues/4033 # moon_separation = Angle([moon.separation(target) for target in targets]) if self.min is None and self.max is not None: mask = self.max > moon_separation elif self.max is None and self.min is not None: mask = self.min < moon_separation elif self.min is not None and self.max is not None: mask = ((self.min < moon_separation) & (moon_separation < self.max)) else: raise ValueError("No max and/or min specified in " "MoonSeparationConstraint.") return mask
[docs]class MoonIlluminationConstraint(Constraint): """ Constrain the fractional illumination of the Earth's moon. """ def __init__(self, min=None, max=None): """ Parameters ---------- min : float or `None` (optional) Minimum acceptable fractional illumination. `None` indicates no limit. max : float or `None` (optional) Maximum acceptable fractional illumination. `None` indicates no limit. """ self.min = min self.max = max
[docs] def compute_constraint(self, times, observer, targets): illumination = np.array(moon_illumination(times, observer.location)) if self.min is None and self.max is not None: mask = self.max > illumination elif self.max is None and self.min is not None: mask = self.min < illumination elif self.min is not None and self.max is not None: mask = ((self.min < illumination) & (illumination < self.max)) else: raise ValueError("No max and/or min specified in " "MoonSeparationConstraint.") return mask
[docs]class LocalTimeConstraint(Constraint): """ Constrain the observable hours. """ def __init__(self, min=None, max=None): """ Parameters ---------- min : `~datetime.time` Earliest local time. `None` indicates no limit. max : `~datetime.time` Latest local time. `None` indicates no limit. Examples -------- Constrain the observations to targets that are observable between 23:50 and 04:08 local time: >>> from astroplan import Observer >>> from astroplan.constraints import LocalTimeConstraint >>> import datetime as dt >>> subaru = Observer.at_site("Subaru", timezone="US/Hawaii") >>> constraint = LocalTimeConstraint(min=dt.time(23,50), max=dt.time(4,8)) # bound times between 23:50 and 04:08 local Hawaiian time """ self.min = min self.max = max if self.min is None and self.max is None: raise ValueError("You must at least supply either a minimum or a maximum time.") if self.min is not None: if not isinstance(self.min, datetime.time): raise TypeError("Time limits must be specified as datetime.time objects.") if self.max is not None: if not isinstance(self.max, datetime.time): raise TypeError("Time limits must be specified as datetime.time objects.")
[docs] def compute_constraint(self, times, observer, targets): timezone = None # get timezone from time objects, or from observer if self.min is not None: timezone = self.min.tzinfo elif self.max is not None: timezone = self.max.tzinfo if timezone is None: timezone = observer.timezone if self.min is not None: min_time = self.min else: min_time = self.min = datetime.time(0, 0, 0) if self.max is not None: max_time = self.max else: max_time = datetime.time(23, 59, 59) # If time limits occur on same day: if self.min < self.max: mask = [min_time < t.datetime.time() < max_time for t in times] # If time boundaries straddle midnight: else: mask = [(t.datetime.time() > min_time) or (t.datetime.time() < max_time) for t in times] return mask
[docs]def is_always_observable(constraints, observer, targets, times=None, time_range=None, time_grid_resolution=0.5*u.hour): """ Are the ``targets`` always observable throughout ``time_range`` given constraints in ``constraints_list`` for ``observer``? Parameters ---------- constraints : list or `~astroplan.constraints.Constraint` Observational constraint(s) observer : `~astroplan.Observer` The observer who has constraints ``constraints`` targets : {list, `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`} Target or list of targets times : `~astropy.time.Time` (optional) Array of times on which to test the constraint time_range : `~astropy.time.Time` (optional) Lower and upper bounds on time sequence, with spacing ``time_resolution``. This will be passed as the first argument into `~astroplan.time_grid_from_range`. time_resolution : `~astropy.units.Quantity` (optional) If ``time_range`` is specified, determine whether constraints are met between test times in ``time_range`` by checking constraint at linearly-spaced times separated by ``time_resolution``. Default is 0.5 hours. Returns ------- ever_observable : list List of booleans of same length as ``targets`` for whether or not each target is observable in the time range given the constraints. """ if not hasattr(constraints, '__len__'): constraints = [constraints] applied_constraints = [constraint(observer, targets, times=times, time_range=time_range, time_grid_resolution=time_grid_resolution) for constraint in constraints] contraint_arr = np.logical_and.reduce(applied_constraints) return np.all(contraint_arr, axis=1)
[docs]def is_observable(constraints, observer, targets, times=None, time_range=None, time_grid_resolution=0.5*u.hour): """ Are the ``targets`` observable during ``time_range`` given constraints in ``constraints_list`` for ``observer``? Parameters ---------- constraints : list or `~astroplan.constraints.Constraint` Observational constraint(s) observer : `~astroplan.Observer` The observer who has constraints ``constraints`` targets : {list, `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`} Target or list of targets times : `~astropy.time.Time` (optional) Array of times on which to test the constraint time_range : `~astropy.time.Time` (optional) Lower and upper bounds on time sequence, with spacing ``time_resolution``. This will be passed as the first argument into `~astroplan.time_grid_from_range`. time_resolution : `~astropy.units.Quantity` (optional) If ``time_range`` is specified, determine whether constraints are met between test times in ``time_range`` by checking constraint at linearly-spaced times separated by ``time_resolution``. Default is 0.5 hours. Returns ------- ever_observable : list List of booleans of same length as ``targets`` for whether or not each target is ever observable in the time range given the constraints. """ if not hasattr(constraints, '__len__'): constraints = [constraints] applied_constraints = [constraint(observer, targets, times=times, time_range=time_range, time_grid_resolution=time_grid_resolution) for constraint in constraints] contraint_arr = np.logical_and.reduce(applied_constraints) return np.any(contraint_arr, axis=1)
[docs]def observability_table(constraints, observer, targets, times=None, time_range=None, time_grid_resolution=0.5*u.hour): """ Creates a table with information about observablity for all the ``targets`` over the requeisted ``time_range``, given the constraints in ``constraints_list`` for ``observer``. Parameters ---------- constraints : list or `~astroplan.constraints.Constraint` Observational constraint(s) observer : `~astroplan.Observer` The observer who has constraints ``constraints`` targets : {list, `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`} Target or list of targets times : `~astropy.time.Time` (optional) Array of times on which to test the constraint time_range : `~astropy.time.Time` (optional) Lower and upper bounds on time sequence, with spacing ``time_resolution``. This will be passed as the first argument into `~astroplan.time_grid_from_range`. time_resolution : `~astropy.units.Quantity` (optional) If ``time_range`` is specified, determine whether constraints are met between test times in ``time_range`` by checking constraint at linearly-spaced times separated by ``time_resolution``. Default is 0.5 hours. Returns ------- observability_table : `~astropy.table.Table` A Table containing the observability information for each of the ``targets``. The table contains four columns with information about the target and it's observability: ``'target name'``, ``'ever observable'``, ``'always observable'``, and ``'fraction of time observable'``. It also contains metadata entries ``'times'`` (with an array of all the times), ``'observer'`` (the `~astroplan.Observer` object), and ``'constraints'`` (containing the supplied ``constraints``). """ if not hasattr(constraints, '__len__'): constraints = [constraints] applied_constraints = [constraint(observer, targets, times=times, time_range=time_range, time_grid_resolution=time_grid_resolution) for constraint in constraints] contraint_arr = np.logical_and.reduce(applied_constraints) colnames = ['target name', 'ever observable', 'always observable', 'fraction of time observable'] target_names = [target.name for target in targets] ever_obs = np.any(contraint_arr, axis=1) always_obs = np.all(contraint_arr, axis=1) frac_obs = np.sum(contraint_arr, axis=1) / contraint_arr.shape[1] tab = table.Table(names=colnames, data=[target_names, ever_obs, always_obs, frac_obs]) if times is None and time_range is not None: times = time_grid_from_range(time_range, time_resolution=time_grid_resolution) tab.meta['times'] = times.datetime tab.meta['observer'] = observer tab.meta['constraints'] = constraints return tab