Lecturer: Robert Quimby
Jupyter Notebook Author: Shubham Srivastav & Cameron Hummels
This is a Jupyter notebook lesson taken from the GROWTH Winter School 2018. 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
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-04-02 05:53:06.289935 2458575.7452116893 58575.24521168906 2019.24998688134
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-04-02 07:38:06.289935 UTC
Update the bulletin:
download_IERS_A()
print("Available observatories: \n{0}"
.format(', '.join(EarthLocation.get_site_names())))
Available observatories: , , , ALMA, Anglo-Australian Observatory, Apache Point, Apache Point Observatory, Atacama Large Millimeter Array, BAO, Beijing XingLong Observatory, Black Moshannon Observatory, CHARA, Canada-France-Hawaii Telescope, Catalina Observatory, Cerro Pachon, Cerro Paranal, Cerro Tololo, Cerro Tololo Interamerican Observatory, DCT, 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, National Observatory of Venezuela, Noto, Observatorio Astronomico Nacional, San Pedro Martir, Observatorio Astronomico Nacional, Tonantzintla, Palomar, Paranal Observatory, Roque de los Muchachos, SAAO, SALT, SRT, Siding Spring Observatory, Southern African Large Telescope, Subaru, Subaru Telescope, Sutherland, TUG, UKIRT, United Kingdom Infrared Telescope, Vainu Bappu Observatory, Very Large Array, W. M. Keck Observatory, Whipple, Whipple Observatory, aao, alma, apo, bmo, cfht, ctio, dao, dct, 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, srt, sso, tona, tug, ukirt, vbo, vla
#Indian Astronomical Observatory is not listed in the database, so let's define the location
longitude = '78d57m53s'
latitude = '32d46m44s'
elevation = 4500 * u.m
location = EarthLocation.from_geodetic(longitude, latitude, elevation)
iaohanle = Observer(location = location, timezone = 'Asia/Kolkata',
name = "IAO", description = "GROWTH-India 70cm telescope")
iaohanle
<Observer: name='IAO', location (lon, lat, el)=(78.96472222222222 deg, 32.77888888888889 deg, 4499.999999999798 m), timezone=<DstTzInfo 'Asia/Kolkata' LMT+5:53:00 STD>>
#Calculating the sunset, midnight and sunrise times for our observatory
#What is astronomical twilight?
sunset_iao = iaohanle.sun_set_time(now, which='nearest')
eve_twil_iao = iaohanle.twilight_evening_astronomical(now, which='nearest')
midnight_iao = iaohanle.midnight(now, which='next')
morn_twil_iao = iaohanle.twilight_morning_astronomical(now, which='next')
sunrise_iao = iaohanle.sun_rise_time(now, which='next')
print("Sunset at IAO will be at {0.iso} UTC".format(sunset_iao))
print("Astronomical evening twilight at IAO will be at {0.iso} UTC".format(eve_twil_iao))
print("Midnight at IAO will be at {0.iso} UTC".format(midnight_iao))
print("Astronomical morning twilight at IAO will be at {0.iso} UTC".format(morn_twil_iao))
print("Sunrise at IAO will be at {0.iso} UTC".format(sunrise_iao))
Sunset at IAO will be at 2019-04-02 13:00:32.540 UTC Astronomical evening twilight at IAO will be at 2019-04-02 14:28:19.957 UTC Midnight at IAO will be at 2019-04-02 18:47:32.580 UTC Astronomical morning twilight at IAO will be at 2019-04-02 23:06:42.282 UTC Sunrise at IAO will be at 2019-04-03 00:34:25.511 UTC
Find the effective length of time (in hours) available for astronomical observations at IAO tonight
observing_time = (morn_twil_iao - eve_twil_iao).to(u.h)
print("You can observe for {0:.1f} at IAO tonight".format(observing_time))
You can observe for 8.6 h at IAO tonight
#What is the LST now at IAO Hanle?
#What would the LST be at IAO at local midnight?
lst_now = iaohanle.local_sidereal_time(now)
lst_mid = iaohanle.local_sidereal_time(midnight_iao)
print("LST at IAO now is {0:.2f}".format(lst_now))
print("LST at IAO at local midnight will be {0:.2f}".format(lst_mid))
LST at IAO now is 23.84 hourangle LST at IAO at local midnight will be 12.78 hourangle
Targets can be defined by name or coordinates.
coords = SkyCoord('00h42m41.8s', '+40d51m55.0s', frame='icrs') # coordinates of Andromeda Galaxy (M32)
m32 = FixedTarget(name = 'M32', coord=coords)
m32.ra.hms
hms_tuple(h=0.0, m=42.0, s=41.799999999999926)
#Check to see if target is up at evening twilight.
#Also check if target is available at midnight and morning twilight.
print(iaohanle.target_is_up(eve_twil_iao, m32))
print(iaohanle.target_is_up(midnight_iao, m32))
print(iaohanle.target_is_up(morn_twil_iao, m32))
True False True
#Altitude and Azimuth of target
aa = iaohanle.altaz(eve_twil_iao, m32)
aa.alt.degree, aa.az.degree
(4.53096054709229, 316.98764800592306)
Checking rise times of targets
m32rise = iaohanle.target_rise_time(now, m32, which = 'next', horizon = 0 * u.deg)
print(m32rise.iso) #default format is JD
2019-04-02 22:27:58.640
target = FixedTarget.from_name('m51') #Messier 51
target.coord
<SkyCoord (ICRS): (ra, dec) in deg (202.469575, 47.1952583)>
get_body('jupiter', now)
<SkyCoord (GCRS: obstime=2019-04-02 05:53:06.289935, obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s): (ra, dec, distance) in (deg, deg, AU) (263.46769437, -22.66809942, 4.92780141)>
#get moon position at midnight
get_moon(midnight_iao)
<SkyCoord (GCRS: obstime=2458576.2830159706, obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s): (ra, dec, distance) in (deg, deg, km) (346.38883828, -10.02541931, 404098.19517621)>
#How bright is the moon at midnight?
moon_illumination(midnight_iao)
0.0635173895581661
#We can turn solar system objects into 'pseudo-fixed' targets to plan observations
saturn_midnight = FixedTarget(name = 'Saturn', coord = get_body('saturn', midnight_iao))
saturn_midnight.coord
<SkyCoord (GCRS: obstime=2458576.2830159706, obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s): (ra, dec, distance) in (deg, deg, AU) (291.20444303, -21.59108981, 10.13197767)>
#Is the target up at IAO at midnight?
iaohanle.target_is_up(midnight_iao, target)
True
#lets check the alt and az of the target at midnight
target_altaz = iaohanle.altaz(midnight_iao, target)
target_altaz.altaz
<SkyCoord (AltAz: obstime=2458576.2830159706, location=(1028191.03516148, 5272251.01802928, 3435803.52567405) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg (27.08622303, 73.41928233)>
That's a good enough elevation to observe the target.
#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, iaohanle, midnight_iao)
/home/chummels/miniconda3/lib/python3.7/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: The frac parameter was deprecated in version 2.1. Use tick padding via Axes.tick_params instead. warnings.warn(message, mplDeprecation, stacklevel=1)
<matplotlib.axes._subplots.PolarAxesSubplot at 0x7f85856100f0>
Now let us see how the target moves over the course of the night
t_start = eve_twil_iao
t_end = morn_twil_iao
t_observe = t_start + (t_end - t_start) * np.linspace(0.0, 1.0, 20)
plot_sky(target, iaohanle, t_observe)
<matplotlib.axes._subplots.PolarAxesSubplot at 0x7f85869f8828>
Now let's plot the airmass as a function of time
plot_airmass(target, iaohanle, t_observe)
plt.ylim(4,0.5)
(4, 0.5)
The airmass is above 2 for the better part of the night, making M51 a good target to observe from IAO tonight. 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
plot_finder_image(target, fov_radius = 20 * u.arcmin) #field of view corresponding to the GROWTH-India telesocpe
(<matplotlib.axes._subplots.WCSAxesSubplot at 0x7f85878706d8>, <astropy.io.fits.hdu.image.PrimaryHDU at 0x7f8588c6e400>)
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.633083, 22.0145)>, <FixedTarget "m42" at SkyCoord (ICRS): (ra, dec) in deg (83.82208, -5.39111)>, <FixedTarget "m55" at SkyCoord (ICRS): (ra, dec) in deg (294.998792, -30.96475)>]
Which of these targets is up now?
iaohanle.target_is_up(now, targets)
array([ True, True, True, True, True])
iaohanle.target_is_up(midnight_iao, targets)
array([ True, True, False, False, False])
Find out the times at which the targets rise to an elevation of 10 degrees. Use target_rise_time.
for target in targets:
print(iaohanle.target_rise_time(now, target, which = 'next', horizon = 10*u.deg).iso)
2019-04-02 17:40:03.863 -4715-02-28 12:00:00.000 2019-04-03 05:26:07.416 2019-04-02 06:41:16.114 2019-04-02 22:16:38.379
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]
What is the elevation of Vega now?
iaohanle.altaz(now, targets[0])
<SkyCoord (AltAz: obstime=2019-04-02 05:53:06.289935, location=(1028191.03516148, 5272251.01802928, 3435803.52567405) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg (299.98263436, 28.28822783)>
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 = iaohanle.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.show()
Plot the altitude as a function of time for tonight for each of the targets in a single plot
times = t_start + (t_end - t_start) * np.linspace(0.0, 1.0, 20)
elevations = []
for target in targets:
elevations.append(iaohanle.altaz(times, target).alt)
ax = plt.gca()
for elevation in elevations:
ax.plot_date(times.plot_date, elevation)
ax.set(xlabel = 'Time UTC [MM-DD HH]' ,ylabel = 'Altitude [deg]')
plt.setp(ax.get_xticklabels(), rotation=45, ha='right')
plt.legend()
No handles with labels found to put in legend.
<matplotlib.legend.Legend at 0x7f85874de160>
Plot sky positions for each target using plot_sky for tonight at IAO 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, iaohanle, times)
plt.legend(loc=[1.1,0])
/home/chummels/miniconda3/lib/python3.7/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: The frac parameter was deprecated in version 2.1. Use tick padding via Axes.tick_params instead. warnings.warn(message, mplDeprecation, stacklevel=1)
<matplotlib.legend.Legend at 0x7f8587488b70>
Plot airmass vs time for each target in targets for tonight at IAO.
for target in targets:
plot_airmass(target, iaohanle, t_observe)
plt.ylim(4,0.5)
plt.legend()
<matplotlib.legend.Legend at 0x7f85873c7438>
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, iaohanle, targets, time_range=t_range)
print(ever_observable)
# Are targets always observable in the time range?
always_observable = is_always_observable(constraints, iaohanle, targets, time_range=t_range)
print(always_observable)
# During what months are the targets ever observable?
obs_months = months_observable(constraints, iaohanle, targets)
[ True True True True False] [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 True False m55 False False
Or we could do this directly using the observability_table function
from astroplan import observability_table
table = observability_table(constraints, iaohanle, targets, time_range = t_range)
print(table)
target name ever observable always observable fraction of time observable ----------- --------------- ----------------- --------------------------- vega True False 0.5 polaris True True 1.0 m1 True False 0.3 m42 True False 0.2 m55 False False 0.0
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