visible_ras.py

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

from astroplan import (Observer, FixedTarget, observability_table,
                       AltitudeConstraint, AtNightConstraint)
from astropy.time import Time
import astropy.units as u
import numpy as np
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt

RAs = np.arange(0, 360+15, 15)*u.deg
Decs = -31*np.ones(len(RAs))*u.deg

c = SkyCoord(ra=RAs, dec=Decs, frame='icrs')
ft = [FixedTarget(coord=c_i) for c_i in c]

constraints = [AltitudeConstraint(min=30*u.deg),
               AtNightConstraint.twilight_nautical()]

obs = Observer.at_site("AAO")

test_dates = Time('2000-01-01') + np.linspace(0, 2*365, 20)*u.day

fraction_matrix = np.zeros((len(RAs), len(test_dates)))

for i, test_date in enumerate(test_dates):
    time_range = test_date + u.Quantity([0, 24])*u.hour

    table = observability_table(constraints, obs, ft, time_range=time_range)
    fraction_matrix[:, i] = table['fraction of time observable'].data

fig, ax = plt.subplots(1, 1, figsize=(10, 10))

ax.imshow(fraction_matrix, interpolation='nearest', origin='lower')
ax.set_xticks((test_dates.unix - test_dates[0].unix)/86400)
ax.set_xticklabels([dt.strftime('%b') for dt in test_dates.datetime])

ax.set_yticks(RAs[::2].to(u.hourangle).value)
ax.set_yticklabels(RAs[::2].to(u.hourangle).value)

ax.set(ylabel='RA')
plt.show()