Lecturer: Robert Quimby
Jupyter Notebook Author: Shubham Srivastav, Cameron Hummels & Robert Quimby
This is a Jupyter notebook lesson taken from the GROWTH Summer School 2019. For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-astro-school-2018-resources.html
Demonstrate how to plan observations prior to an observing run.
See GROWTH school webpage for detailed instructions on how to install these modules and packages. Nominally, you should be able to install the python modules with pip install <module>
. The external astromatic packages are easiest installed using package managers (e.g., rpm
, apt-get
).
None
import numpy as np
from astropy import units as u
from astropy.time import Time
from astropy.coordinates import SkyCoord
from astropy.coordinates import EarthLocation
import pytz
%matplotlib inline
from astroplan import Observer, FixedTarget
from astropy.utils.iers import conf
conf.auto_max_age = None
from astroplan import download_IERS_A
from astropy.coordinates import get_sun, get_moon, get_body
from astroplan import moon_illumination
WARNING: AstropyDeprecationWarning: astropy.extern.six will be removed in 4.0, use the six module directly if it is still needed [astropy.extern.six]
date = Time("2018-12-03", format='iso')
print(date)
2018-12-03 00:00:00.000
now = Time.now()
print(now)
print(now.jd)
print(now.mjd)
print(now.decimalyear)
2019-07-09 09:28:32.298729 2458673.8948182724 58673.394818272325 2019.5188899130749
What time will it be (in UTC) after 1 hour 45 minutes from now
? Complete the line below to print it out.
print("In 1 hour and 45 minutes, the time will be {0} UTC".format(now + 1*u.h + 45*u.min))
In 1 hour and 45 minutes, the time will be 2019-07-09 11:13:32.298729 UTC
Update the bulletin:
download_IERS_A()
print("Available observatories: \n{0}"
.format(', '.join(EarthLocation.get_site_names())))
Available observatories: , , , ALMA, ATST, Anglo-Australian Observatory, Apache Point, Apache Point Observatory, Atacama Large Millimeter Array, BAO, BBSO, Beijing XingLong Observatory, Black Moshannon Observatory, CHARA, Canada-France-Hawaii Telescope, Catalina Observatory, Cerro Pachon, Cerro Paranal, Cerro Tololo, Cerro Tololo Interamerican Observatory, DCT, DKIST, Discovery Channel Telescope, Dominion Astrophysical Observatory, GBT, Gemini South, Green Bank Telescope, Hale Telescope, Haleakala Observatories, Happy Jack, IAO, JCMT, James Clerk Maxwell Telescope, Jansky Very Large Array, Keck Observatory, Kitt Peak, Kitt Peak National Observatory, La Silla Observatory, Large Binocular Telescope, Las Campanas Observatory, Lick Observatory, Lowell Observatory, MWA, Manastash Ridge Observatory, McDonald Observatory, Medicina, Medicina Dish, Michigan-Dartmouth-MIT Observatory, Mount Graham International Observatory, Mt Graham, Mt. Ekar 182 cm. Telescope, Mt. Stromlo Observatory, Multiple Mirror Telescope, Murchison Widefield Array, NOV, NST, National Observatory of Venezuela, Noto, Observatorio Astronomico Nacional, San Pedro Martir, Observatorio Astronomico Nacional, Tonantzintla, Palomar, Paranal Observatory, Roque de los Muchachos, SAAO, SALT, SPO, SRT, Sac Peak, Sacramento Peak, Siding Spring Observatory, Southern African Large Telescope, Subaru, Subaru Telescope, Sunspot, Sutherland, TUG, UKIRT, United Kingdom Infrared Telescope, Vainu Bappu Observatory, Very Large Array, W. M. Keck Observatory, Whipple, Whipple Observatory, aao, alma, apo, bbso, bmo, cfht, ctio, dao, dct, dkist, ekar, example_site, flwo, gbt, gemini_north, gemini_south, gemn, gems, greenwich, haleakala, iao, irtf, jcmt, keck, kpno, lapalma, lasilla, lbt, lco, lick, lowell, mcdonald, mdm, medicina, mmt, mro, mso, mtbigelow, mwa, mwo, noto, ohp, paranal, salt, sirene, spm, spo, srt, sso, tona, tug, ukirt, vbo, vla
# Mount Laguna Observatory is not listed in the database, so let's define the location
latitude = 32.842167 * u.deg
longitude = -116.426938 * u.deg
elevation = 1860 * u.m
location = EarthLocation.from_geodetic(longitude, latitude, elevation)
mlo = Observer(location = location, timezone = 'America/Los_Angeles',
name = "MLO", description = "MLO 1.0-m telescope")
mlo
<Observer: name='MLO', location (lon, lat, el)=(-116.426938 deg, 32.842166999999996 deg, 1860.0000000008795 m), timezone=<DstTzInfo 'America/Los_Angeles' LMT-1 day, 16:07:00 STD>>
# ##### just for testing...REMOVE! #####
now = Time('2019-08-04 21:17:59')
# Calculating the sunset, midnight and sunrise times for our observatory
# What is astronomical twilight?
sunset_mlo = mlo.sun_set_time(now, which='nearest')
eve_twil_mlo = mlo.twilight_evening_astronomical(now, which='nearest')
midnight_mlo = mlo.midnight(now, which='next')
morn_twil_mlo = mlo.twilight_morning_astronomical(now, which='next')
sunrise_mlo = mlo.sun_rise_time(now, which='next')
print("Sunset at MLO will be at {0.iso} UTC".format(sunset_mlo))
print("Astronomical evening twilight at MLO will be at {0.iso} UTC".format(eve_twil_mlo))
print("Midnight at MLO will be at {0.iso} UTC".format(midnight_mlo))
print("Astronomical morning twilight at MLO will be at {0.iso} UTC".format(morn_twil_mlo))
print("Sunrise at MLO will be at {0.iso} UTC".format(sunrise_mlo))
Sunset at MLO will be at 2019-08-05 02:37:31.424 UTC Astronomical evening twilight at MLO will be at 2019-08-05 04:14:21.609 UTC Midnight at MLO will be at 2019-08-05 07:51:56.122 UTC Astronomical morning twilight at MLO will be at 2019-08-05 11:29:30.653 UTC Sunrise at MLO will be at 2019-08-05 13:06:21.359 UTC
Find the effective length of time (in hours) available for optical astronomical observations at MLO tonight
observing_time = (morn_twil_mlo - eve_twil_mlo).to(u.hour)
print("You can observe for {0:.1f} at MLO tonight".format(observing_time))
You can observe for 7.3 h at MLO tonight
#What is the LST now at MLO?
#What would the LST be at MLO at local midnight?
lst_now = mlo.local_sidereal_time(now)
lst_mid = mlo.local_sidereal_time(midnight_mlo)
print("LST at MLO now is {0:.2f}".format(lst_now))
print("LST at MLO at local midnight will be {0:.2f}".format(lst_mid))
LST at MLO now is 10.41 hourangle LST at MLO at local midnight will be 21.01 hourangle
Targets can be defined by name or coordinates.
# using coordinates
coords = SkyCoord('18h53m35.097s +33d01m44.8831s', frame='icrs') # coordinates of the Ring Nebula (M57)
m57 = FixedTarget(name = 'M57', coord=coords)
m57.ra.hms
hms_tuple(h=18.0, m=53.0, s=35.09699999999441)
# by name
target = FixedTarget.from_name('m57') # Messier 57
target.coord
<SkyCoord (ICRS): (ra, dec) in deg (283.39623732, 33.02913421)>
Check to see if target is "up" at evening twilight (assume "up" means more than 30 degrees above the horizon). Also check if target is available at midnight and morning twilight.
# check if the target is up
print(mlo.target_is_up(eve_twil_mlo, m57, horizon=30*u.deg))
print(mlo.target_is_up(midnight_mlo, m57, horizon=30*u.deg))
print(mlo.target_is_up(morn_twil_mlo, m57, horizon=30*u.deg))
True True False
# Altitude and Azimuth of target at evening twilight
aa = mlo.altaz(eve_twil_mlo, m57)
aa.alt.degree, aa.az.degree
(70.74681163016001, 83.05825913756574)
Determine the time at which the target rises
m57rise = mlo.target_rise_time(now, m57, which = 'next', horizon=0*u.deg)
print(m57rise.iso) #default format is JD
2019-08-04 22:07:56.216
get_body('jupiter', now)
<SkyCoord (GCRS: obstime=2019-08-04 21:17:59.000, obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s): (ra, dec, distance) in (deg, deg, AU) (253.0253309, -22.09201317, 4.66743262)>
# get moon position at midnight
get_moon(midnight_mlo)
<SkyCoord (GCRS: obstime=2458700.827732899, obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s): (ra, dec, distance) in (deg, deg, km) (192.07868451, 0.46596914, 366500.67856505)>
# How bright is the moon at midnight?
moon_illumination(midnight_mlo)
0.24118319358448292
# We can turn solar system objects into 'pseudo-fixed' targets to plan observations
saturn_midnight = FixedTarget(name = 'Saturn', coord = get_body('saturn', midnight_mlo))
saturn_midnight.coord
<SkyCoord (GCRS: obstime=2458700.827732899, obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s): (ra, dec, distance) in (deg, deg, AU) (286.35261842, -22.30625221, 9.13525905)>
#Is the target up at MLO at midnight?
mlo.target_is_up(midnight_mlo, target)
True
#lets check the alt and az of the target at midnight
target_altaz = mlo.altaz(midnight_mlo, target)
target_altaz.altaz
<SkyCoord (AltAz: obstime=2458700.827732899, location=(-2387971.90799836, -4804864.68705021, 3440274.03123442) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg (279.19059776, 63.62151349)>
That is high over head; ideal for observing.
#Find the airmass
target_altaz.secz
Now we can visualize what we have done so far using some plots
import matplotlib.pyplot as plt
from astroplan.plots import plot_sky, plot_airmass
#position of target at midnight
plot_sky(target, mlo, midnight_mlo);
Now let us see how the target moves over the course of the night
t_start = eve_twil_mlo
t_end = morn_twil_mlo
t_observe = t_start + (t_end - t_start) * np.linspace(0.0, 1.0, 20)
plot_sky(target, mlo, t_observe);
Now let's plot the airmass as a function of time
plot_airmass(target, mlo, t_observe)
plt.grid();
The airmass is above 2 for the better part of the night, making M57 a good summer target from MLO. Note that the default airmass limit is 3 in astroplan, corresponding to ~19 degrees elevation.
from astroplan.plots import plot_finder_image
from astroquery.skyview import SkyView
Load an image of the field in which the target lies.
# field of view corresponding to the MLO 1.0-m telesocpe
fov = 14*u.arcmin
# plot the image
plot_finder_image(target, fov_radius=fov);
Now let's define an array of targets to work with
target_names = ['vega', 'polaris', 'm1', 'm42', 'm55']
targets = [FixedTarget.from_name(target) for target in target_names]
targets
[<FixedTarget "vega" at SkyCoord (ICRS): (ra, dec) in deg (279.23473479, 38.78368896)>, <FixedTarget "polaris" at SkyCoord (ICRS): (ra, dec) in deg (37.95456067, 89.26410897)>, <FixedTarget "m1" at SkyCoord (ICRS): (ra, dec) in deg (83.63308333, 22.0145)>, <FixedTarget "m42" at SkyCoord (ICRS): (ra, dec) in deg (83.82208333, -5.39111111)>, <FixedTarget "m55" at SkyCoord (ICRS): (ra, dec) in deg (294.99879167, -30.96475)>]
Which of these targets are up now?
mlo.target_is_up(now, targets)
array([False, True, True, True, False])
Which of these targets are up at local midnight?
mlo.target_is_up(midnight_mlo, targets)
array([ True, True, False, False, True])
Find out the times at which the targets rise to an elevation of 10 degrees. Use target_rise_time.
for target in targets:
print(mlo.target_rise_time(now, target, which = 'next', horizon = 10*u.deg).iso)
2019-08-04 22:31:41.707 -4715-02-28 12:00:00.000 2019-08-05 10:17:54.688 2019-08-05 11:29:17.468
WARNING: TargetAlwaysUpWarning: Target with index 0 does not cross horizon=10.0 deg within 24 hours [astroplan.observer] WARNING: ErfaWarning: ERFA function "d2dtf" yielded 1 of "dubious year (Note 5)" [astropy._erfa.core]
2019-08-05 03:08:54.426
How high is Vega above the horizion now?
mlo.altaz(now, targets[0])
<SkyCoord (AltAz: obstime=2019-08-04 21:17:59.000, location=(-2387971.90799836, -4804864.68705021, 3440274.03123442) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg (40.71340897, -1.06248615)>
Now let's plot the elevation of Vega to see how it varies over the night
times = (t_start - 0.5 * u.h) + (t_end - t_start + 1 * u.h) * np.linspace(0.0, 1.0, 40)
elevations = mlo.altaz(times, targets[0]).alt
ax = plt.gca()
ax.plot_date(times.plot_date, elevations.deg)
ax.set(xlabel = 'Time UTC [MM-DD HH]', ylabel = 'Altitude [deg]')
plt.setp(ax.get_xticklabels(), rotation=45, ha='right')
plt.grid()
Plot the altitude as a function of time for tonight for each of the targets in a single plot
ax = plt.gca()
times = t_start + (t_end - t_start) * np.linspace(0.0, 1.0, 20)
for target in targets:
elevation = mlo.altaz(times, target).alt
ax.plot_date(times.plot_date, elevation, label=target.name)
ax.set(xlabel = 'Time UTC [MM-DD HH]', ylabel = 'Altitude [deg]')
plt.setp(ax.get_xticklabels(), rotation=45, ha='right')
plt.legend()
plt.grid()
Plot sky positions for each target using plot_sky for tonight at MLO in a single plot.
times = (t_start - 0.5 * u.h) + (t_end - t_start + 1 * u.h) * np.linspace(0.0, 1.0, 20)
for target in targets:
plot_sky(target, mlo, times)
plt.legend(loc=[1.1,0]);
Plot airmass vs time for each target in targets for tonight at MLO.
for target in targets:
plot_airmass(target, mlo, t_observe)
plt.ylim(4,0.5)
plt.legend()
plt.grid()
You can set specific constraints that define when a target is "observable"
from astroplan import (AltitudeConstraint, AirmassConstraint,
AtNightConstraint, MoonSeparationConstraint)
constraints = [AltitudeConstraint(15*u.deg, 84*u.deg),
AirmassConstraint(3), AtNightConstraint.twilight_civil(), MoonSeparationConstraint(min = 10 * u.deg)]
t_range = Time([t_start - 0.5 * u.hour, t_end + 0.5 * u.hour])
from astroplan import is_observable, is_always_observable, months_observable
# Are targets ever observable in the time range?
ever_observable = is_observable(constraints, mlo, targets, time_range=t_range)
print(ever_observable)
# Are targets always observable in the time range?
always_observable = is_always_observable(constraints, mlo, targets, time_range=t_range)
print(always_observable)
[ True True True False True] [False True False False False]
The functions is_observable and ever_observable return boolean arrays. Let's print their output in tabular form.
from astropy.table import Table
observability_table = Table()
observability_table['targets'] = [target.name for target in targets]
observability_table['ever_observable'] = ever_observable
observability_table['always_observable'] = always_observable
print(observability_table)
targets ever_observable always_observable ------- --------------- ----------------- vega True False polaris True True m1 True False m42 False False m55 True False
Or we could do this directly using the observability_table function
from astroplan import observability_table
table = observability_table(constraints, mlo, targets, time_range = t_range)
print(table)
target name ever observable always observable fraction of time observable ----------- --------------- ----------------- --------------------------- vega True False 0.9411764705882353 polaris True True 1.0 m1 True False 0.11764705882352941 m42 False False 0.0 m55 True False 0.47058823529411764
# During what months are the targets ever observable?
months_observable(constraints, mlo, targets)
[{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}, {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}, {1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12}, {1, 2, 3, 4, 8, 9, 10, 11, 12}, {3, 4, 5, 6, 7, 8, 9, 10, 11}]
from astropy.io import ascii
table = ascii.read('data/targetlist.txt')
targets = [FixedTarget(coord=SkyCoord(ra=ra*u.deg, dec=dec*u.deg), name=name) for name, ra, dec in table]
#The recipe for the remaining part of the exercise is in the previous solved exercises