astropy / astropy/astroplan

LST rise set bar plots

Open
#277 2 comments 0 reactions 0 assignees View on GitHub
enhancement proposal
Dominant language
Python
Stars
225
Forks
121
Avg merge
28m
Merged PRs (30d)
1

Description

I've put together some code to create LST rise/set bar plots like those common at radio observatories. See attached [example](https://github.com/astropy/astroplan/files/677875/20101112_observability.pdf). 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)

````

Contributor guide

No contributing guide indexed for this repository

Research direction

The issue provides a standalone Python example and imports the existing plot_airmass entry point, but names no target file or test. Start by reviewing that plotting entry point and determining where a general LST rise/set plot belongs. Done means the described rise/set, transit, twilight, and LST-axis behavior is available through a reusable routine.

Written by the indexing model from the issue text.

Assessment

Tech stack
matplotlib, numpy, python
Domain
data-visualization
Issue type
Feature
Difficulty
5/5
Estimated time
Over a week
Activity status
Stale
Clarity
Needs clarification
Newbie friendliness
35/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.