LST rise set bar plots

[aakepley]

I've put together some code to create LST rise/set bar plots like those common at radio observatories. See attached example. I don't have time to make it fully general, but feel free to use as the basis for a more general routine.

from astroplan import Observer, FixedTarget
from astroplan.plots import plot_airmass

from astropy.coordinates import SkyCoord, EarthLocation
from astropy.time import Time
from astropy.table import Table
import astropy.units as u

import matplotlib.pyplot as plt
from matplotlib import dates

import numpy as np
import operator

# setup observing
location = EarthLocation.from_geodetic(lon='-79d50m23.406s',lat='38d25m59.236')
gbt = Observer(location=location, elevation=807.43*u.m,name='gbt',timezone='US/Eastern')
time = gbt.midnight(Time('2017-06-01 5:00:00'))
sunset_tonight = gbt.sun_set_time(time, which='nearest')
sunrise_tonight = gbt.sun_rise_time(time, which='nearest')

# set up targets
target_list = Table.read('gbt_target_coord.txt',format='ascii')

# set up plot
fig = figure(figsize=(8.5,11),edgecolor='white')
ax = plt.gca()
time_ax = time + [-12,12]*u.hour
ax.set_xlim([time_ax[0].plot_date,time_ax[-1].plot_date])
ax.set_ylim([0,len(target_list)+1])
ax.yaxis.set_ticklabels('')
ax.axvline(time.plot_date,color='black',linestyle=':')
ax.set_ylabel("")
ax.set_xlabel("Time centered on {0} {1} [UTC]".format(time.datetime.date(), time.datetime.time()))

# Calculate and order twilights and set plotting alpha for each
start = time_ax[0].datetime
twilights = [
    (gbt.sun_set_time(Time(start), which='next').datetime, 0.0),
    (gbt.twilight_evening_civil(Time(start), which='next').datetime, 0.1),
    (gbt.twilight_evening_nautical(Time(start), which='next').datetime, 0.2),
    (gbt.twilight_evening_astronomical(Time(start), which='next').datetime, 0.3),
    (gbt.twilight_morning_astronomical(Time(start), which='next').datetime, 0.4),
    (gbt.twilight_morning_nautical(Time(start), which='next').datetime, 0.3),
    (gbt.twilight_morning_civil(Time(start), which='next').datetime, 0.2),
    (gbt.sun_rise_time(Time(start), which='next').datetime, 0.1),
    ]

twilights.sort(key=operator.itemgetter(0))
for i, twi in enumerate(twilights[1:], 1):
    ax.axvspan(twilights[i - 1][0], twilights[i][0], ymin=0, ymax=1, color='grey', alpha=twi[1])
    
i=1
for target in target_list:
    coordinates = SkyCoord(target['RA']*u.deg, target['Dec']*u.deg, frame='fk5')
    galaxy = FixedTarget(name=target['Name'], coord=coordinates)

    galaxy_rise = gbt.target_rise_time(time, galaxy,horizon=30*u.deg)
    galaxy_set = gbt.target_set_time(time, galaxy,horizon=30*u.deg)
    galaxy_transit = gbt.target_meridian_transit_time(time,galaxy)
    
    if ( (galaxy_set > time_ax[1]) or
         (galaxy_rise < time_ax[0]) or
         (galaxy_rise > galaxy_set)):
        
        ax.plot_date([galaxy_rise.plot_date,time_ax[1].plot_date],[i,i],linestyle='solid',color='black',lw=2,marker="None")

        ax.plot_date([time_ax[0].plot_date,galaxy_set.plot_date],[i,i],linestyle='solid',color='black',lw=2,marker="None")

    else:
    
        ax.plot_date([galaxy_rise.plot_date,galaxy_set.plot_date],[i,i],linestyle='solid',color='black',lw=2,marker="None")
        ax.plot_date([galaxy_rise.plot_date,galaxy_set.plot_date],[i,i],linestyle='None',color='black',marker="|",markersize=10)

    ax.plot_date([galaxy_transit.plot_date],[i],marker='|',color='black',ms=5.0)
    ax.plot_date([galaxy_rise.plot_date],[i],marker='|',color='black',ms=5.0)
    ax.plot_date([galaxy_set.plot_date],[i],marker='|',color='black',ms=5.0)
    
    ax.annotate(target['Name'],(galaxy_transit.plot_date,i),xytext=(0,5),xycoords='data',textcoords='offset points',horizontalalignment='center')
  
    i= i + 1

# add LST axis -- surprisingly annoying. This is hack. I'm presending that sidereal time is UTC time.
utcticks = Time(dates.num2date(ax.get_xticks()))
lstticks = gbt.local_sidereal_time(utcticks)
utclim = Time(dates.num2date(ax.get_xlim()))
lstlim = gbt.local_sidereal_time(utclim)

start_day = time_ax[0].iso.split(' ')[0]
start_lst = Time(start_day+' '+ lstlim[0].to_string(decimal=False,sep=':',precision=2),format='iso',scale='utc')

end_day = time_ax[1].iso.split(' ')[0]
end_lst = Time(end_day+' '+ lstlim[1].to_string(decimal=False,sep=':',precision=2),format='iso',scale='utc')

ax2 = ax.twiny()
ax2.xaxis_date()
ax2.set_xlim([start_lst.plot_date,end_lst.plot_date])

for day in (start_day,end_day):
    start_exclusion = Time(day+' '+ '14:00:00',format='iso',scale='utc')
    end_exclusion = Time(day+' '+ '20:00:00',format='iso',scale='utc')
    ax2.axvspan(start_exclusion.plot_date,end_exclusion.plot_date,facecolor='red',alpha=0.2) #LST exclusion zone

ax2.set_xlabel('Time centered on {0} [LST]'.format(gbt.local_sidereal_time(time).to_string(decimal=False,sep=':',precision=0)))

ax.figure.canvas.draw()

figlab = time.to_datetime().isoformat() + '.pdf'

fig.savefig(figlab)

[bmorris3]

Hey @aakepley – thanks for this! Would you be able to share a version of the file 'gbt_target_coord.txt' so I can run this example and generalize it?

[aakepley]

Sure. Here it is: gbt_target_coord.txt

Note that sometimes the code above doesn't do the wrapping quite right for some start dates.

Let me know if you need any testing effort.