Interactive KBO Notebook¶
In this notebook, we will make a bare-bones application to interactively explore the movement and light profile of KBOs and other moving objects. To fetch data, we will use the Solar System Object Image Search tool provided by the CADC.
To launch this notebook interactively, click the Binder button.
Warning: Binder instances have limited computational resources. A notebook will shutdown after 10 minutes of inactivity and has a maximum of 2 GB of RAM. Additionally, each repository has a maximum of 100 simultaneous users.
Table of Contents¶
1. Introduction¶
The idea for this notebook is to give an outline of a lightweight, interactive tool in a jupyter notebook. We want the ability to display and analyze data interactively. The end goal is an application that allows a user to select a subset of SSOIS images using a date range slider and plot the subset, then select a circular area on the interactive plot and calculate the aperture sum.
2. Setup¶
This tutorial will go through some of the basic functionalities of the python packages such as:
- astropy for astronomical operations
- ipyaladin for a widget for Aladin Lite
- ipywidgets to build interactive widgets
- matplotlib to plot data
- ipympl for the interactive
%matplotlib widgetbackend
3. SSOIS Query¶
In order to get images of a moving object, we will use the Solar System Object Image Search tool provided by the CADC. We will search by object name between a start and end date. To choose an object to search, we will use the Minor Planet Center search tool and search the database where the minimum semi-major axis value is 50 AU. From the results table, we chose the KBO 2006 WG206 to query for. We will search between the dates January 18th, 2007, to January 19th, 2007.
%matplotlib widget
from urllib.parse import urlencode
import pandas as pd
# Define SSOIS search parameters
kbo_name = '2006 WG206'
start_date = '2007 01 18'
end_date = '2007 01 19'
min_exp_time = 20
# Build query url
params = {
'lang': 'en',
'object': kbo_name,
'search': 'bynameCADC',
'epoch1': start_date,
'epoch2': end_date,
'eunits': 'arcseconds',
'extres': 'no', # Return the extension where the KBO can be seen in
'xyres': 'no', # Return the pixel X and Y position of the KBO
'format': 'tsv'
}
base_url = 'https://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/cadcbin/ssos/ssosclf.pl'
url = '{}?{}'.format(base_url, urlencode(params))
# Access data and convert to panadas dataframe
data_table = pd.read_csv(url, sep='\t')
# Get results with a long exposure and the g-band filter
data_table = data_table[(data_table['Exptime'] > min_exp_time)
& (data_table['Filter'] == "G.MP9401")]
data_table
3.1 Aladin Widget¶
To show the location of the SSOIS results, we can use the interactive widget provided by Aladin. This widget has much more functionality which can be explored using the jupyter notebooks on the ipyaladin github.
import numpy as np
from astropy.table import Column, Table
from ipyaladin import Aladin
from ipywidgets import Layout
# Convert data table to astropy table and rename RA and Dec columns
table = Table.from_pandas(data_table)
table.rename_column('Object_RA', 'RA')
table.rename_column('Object_Dec', 'DEC')
table.rename_column('Telescope/Instrument', 'Telescope_Instrument')
table['Image'] = Column(data=data_table['Image'],
name='Image',
dtype=np.dtype('S10'))
target_pos = '{} {}'.format(table[0]['RA'], table[0]['DEC'])
# Display the widget and add the table
aladin = Aladin(target=target_pos,
layout=Layout(width='100%', height='400px'),
fov=1)
display(aladin)
aladin.add_table(table)
3.2 Get CADC data¶
Since we want to later use the datalink service to fetch cutouts of the KBOs, we will query CADC for the results to get the publisherIDs. In this case, we upload the SSOIS results table to the server and perform a join on the image product ID with the Plane and Observation tables. This will return all the results from the SSOIS query filled with the Plane and Observation data from the CADC which will later be used to fetch the data.
from astroquery.cadc import Cadc
cadc = Cadc()
cadc_results = cadc.exec_sync(
'''SELECT upload_table.*, Plane.publisherID
FROM tap_upload.upload_table AS upload_table
JOIN caom2.Plane AS Plane ON Plane.productID = upload_table.Image
JOIN caom2.Observation AS Observation ON Plane.obsID = Observation.obsID ''',
uploads={'upload_table': table[['Image', 'MJD', 'RA', 'DEC']]})
cadc_table = cadc_results.to_pandas()
# Decode all byte string results to utf-8
str_results = cadc_table[['Image', 'publisherID'
]].stack().str.decode('utf-8').unstack()
for col in str_results:
cadc_table[col] = str_results[col]
cadc_table
4. Helper Functions¶
The end goal is to make an application with a date slider that selects a subset of observations and then displays them, with interactivity to preform aperture photometry. To improve performance, we use image cutouts to display only the chunk of the image that is expected to contain the target of interest. In order to do all of this, we need helper functions to select a subset of images, calculate the cutout parameters, get the cutout urls, and fetch the cutout data.
In the following functions,
get_selection()will be used by the date range slider to return the SSOIS results in the selected date range.get_data_list()will be used the the replot button to calculate the cutout parameters for all the images in the table (using_get_cutout_params()), download the cutouts, then add all the image data to a list of dictionaries (using_get_data_dict())
import math
import re
from urllib.error import HTTPError
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.io import fits
from astropy.time import Time
from astropy.wcs import WCS
from astropy.wcs.utils import skycoord_to_pixel
def get_selection(table, date_range):
"""Returns the subset of rows in the table within the date_range."""
return table[(table['MJD'] >= date_range[0])
& (table['MJD'] < date_range[1])]
def _get_cutout_params(table, min_radius=0.005):
"""Given a table with R.A. and Dec. columns, returns the parameters of the
circle needed to include all positions."""
# Get the max and min R.A. and Dec. values
ra_max, ra_min = table['RA'].max(), table['RA'].min()
dec_max, dec_min = table['DEC'].max(), table['DEC'].min()
# Get the R.A. and Dec. centre and radius that covers all points
ra_centre = (ra_max + ra_min) / 2
dec_centre = (dec_max + dec_min) / 2
ra_radius = math.cos(math.radians(dec_centre)) * abs(ra_max - ra_min) / 2
dec_radius = abs(dec_max - dec_min) / 2
coords = SkyCoord(ra_centre, dec_centre, unit=(u.deg, u.deg))
radius = max(round(max(ra_radius, dec_radius), 2), min_radius) * u.deg
return coords, radius
def _get_data_dict(hdu, ra, dec, date, pix_thres=10):
"""Helper function to return a dict of info for each HDU"""
data = hdu.data
if data is not None:
wcs = WCS(hdu.header)
xp, yp = skycoord_to_pixel(SkyCoord(ra, dec, unit=(u.deg, u.deg)), wcs)
y_max, x_max = data.shape
if xp > pix_thres and x_max - xp > pix_thres:
date_fmtd = Time(date, format='mjd', out_subfmt='date_hm').iso
return {
'image_data': data,
'wcs': wcs,
'ra': ra,
'dec': dec,
'date': date_fmtd
}
def get_data_list(table, ext=0, verbose=False):
"""Given the table with cutout urls, return the fits image and WCS information"""
data_list = []
if len(table) < 1:
print("Range is too small")
return None
coords, radius = _get_cutout_params(table)
table['cutout_url'] = cadc.get_image_list(table, coords, radius)
for url, ra, dec, date in zip(table['cutout_url'], table['RA'],
table['DEC'], table['MJD']):
try:
with fits.open(url, ignore_missing_end=True) as hdulist:
if verbose:
hdulist.info()
for hdu in hdulist:
data_dict = _get_data_dict(hdu, ra, dec, date)
if data_dict:
data_list.append(data_dict)
except HTTPError as ex:
print(ex)
continue
return data_list
5. Widgets¶
Now we make our widgets! We want one widget to allow the user to select a date range to display, a button to replot the data, and another widget to display the aperture data. Technically, we also want the plot itself to be a widget, but we will come to that later.
import math
import numpy as np
from ipywidgets import widgets
from traitlets import traitlets
# Date range slider widget
def make_date_range_slider(table):
""" Function to generate the date range slider based on the date
ranges in the results table"""
hour = 1. / 24.
# Round the start date and end dates for the slider
start_date = hour * math.floor(table['MJD'].min() / hour)
end_date = hour * math.ceil(table['MJD'].max() / hour)
dates = Time(np.arange(start_date, end_date + hour, hour),
format='mjd',
out_subfmt='date_hm')
# Match value to formated date output for nicer display
options = [(' {} '.format(date.iso), date.mjd) for date in dates]
# Range slider to select dates
range_slider = widgets.SelectionRangeSlider(options=options,
index=(0, 1),
description='Date Range:',
disabled=False,
continuous_update=True,
readout=False,
layout={'width': '300px'},
style={
'description_width':
'initial',
'contenteditable': 'false'
})
range_label = widgets.Label(value='{} - {}'.format(*range_slider.label))
# The slider controls the selected rows of the table data
selection = get_selection(table, range_slider.value)
range_slider.add_traits(selection=traitlets.Any(selection))
# When slider changes, update selection table and results num label
def update_selection(change):
new_selection = get_selection(table, change['new'])
range_slider.selection = new_selection
range_label.value = '{} - {}'.format(*range_slider.label)
range_slider.observe(update_selection, names='value')
return widgets.HBox([range_slider, range_label])
# Replot button widget
def make_replot_button(range_slider, kbo_fig, kbo_name):
""" Function to generate the replot button to update the
KBO figure"""
replot_button = widgets.Button(description="Replot")
# When button is clicked, get new selection and replot
def on_button_clicked(change):
data_list = get_data_list(Table.from_pandas(
range_slider.children[0].selection),
ext=-1)
if data_list is not None:
plot_data(kbo_fig, data_list, kbo_name)
replot_button.on_click(on_button_clicked)
return replot_button
# Aperture results widget
def make_aperture_text_widget():
"""Returns widget to display aperture information"""
aperture_widget = widgets.HTML(
value='Click and drag on a subplot to select an area to analyze',
placeholder='<b>Aperture Sum: </b>',
description='<b>Aperture Sum: </b>',
disabled=False,
style={'description_width': 'initial'},
layout=Layout(width='90%'),
)
return aperture_widget
5.1 The KBO Figure¶
This is the interactive plot of the KBO. Since the matplotlib widget backend is used, the figure canvas can be used as a widget, which we will use in a box along with the aperture widget. The figure should plot each object in the data list on axes that share a common WCS. When one axis is moved, the rest should follow. The approximate position of the KBO should be outlined for easier identification. We also want the ability to draw circles to show which area is going to be used in the aperture analysis, meaning we will have to use click event handlers.
import warnings
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
from astropy.visualization import ImageNormalize, LinearStretch, ZScaleInterval
from matplotlib.patches import Circle
from photutils import CircularAperture, SkyAperture, aperture_photometry
# Supress fits processing warnings
warnings.simplefilter('ignore')
def kbo_figure(data_list, kbo_name):
"""Build the figure and click events"""
global data_dict
p_center = None
data_dict = {}
def onclick(event):
global p_center
# Return if in pan/zoom mode or click outside of axes
if (plt.get_current_fig_manager().toolbar.mode != ''
or event.inaxes is None):
return
p_center = (event.xdata, event.ydata)
def onrelease(event):
global p_center
global kbo_fig
global data_dict
# Return if in pan/zoom mode or click outside of axes
if (plt.get_current_fig_manager().toolbar.mode != ''
or event.inaxes is None):
return
# Calculate circle radius
p_outside = (event.xdata, event.ydata)
radius = math.sqrt((p_center[0] - p_outside[0])**2 +
(p_center[1] - p_outside[1])**2)
if radius == 0.0:
return
# Remove old circles
for ax in kbo_fig.axes:
for artist in ax.artists:
if isinstance(artist, plt.Circle):
artist.remove()
# Draw new circles
circle = plt.Circle(p_center, radius, color='black', fill=False)
ax = event.inaxes
ax.add_artist(circle)
kbo_fig.canvas.draw()
# Do aperture photometry
aperture = CircularAperture(p_center, r=radius)
image_data = data_dict[(ax.rowNum, ax.colNum)]
phot_table = aperture_photometry(image_data, aperture)
aperture_widget.value = '{:.8}'.format(phot_table['aperture_sum'][0])
# Create figure and oonnect click events
kbo_fig = plt.figure(figsize=(9, 8))
cid_press = kbo_fig.canvas.mpl_connect('button_press_event', onclick)
cid_release = kbo_fig.canvas.mpl_connect('button_release_event', onrelease)
plot_data(kbo_fig, data_list, kbo_name)
return kbo_fig
def plot_data(kbo_fig, data_list, kbo_name, max_cols=3):
"""Plot data on datalist to KBO fig"""
global data_dict
data_dict = dict()
data_list = sorted(data_list, key=lambda data: data['date'])
n = len(data_list)
# Clear the figure if needed
plt.clf()
kbo_fig.suptitle('Images of KBO {}'.format(kbo_name),
fontsize=14,
fontweight='bold')
# Set all axes to be in same WCS
wcs_trans = data_list[0]['wcs']
gs = gridspec.GridSpec(math.ceil(n / max_cols),
max_cols,
figure=kbo_fig,
wspace=0.1)
for idx, data in enumerate(data_list):
try:
ax = plt.subplot(gs[idx // max_cols, idx % max_cols],
projection=wcs_trans,
adjustable='box',
aspect='equal')
# Set the global data_dict to be used for photometry
data_dict[(ax.rowNum, ax.colNum)] = data['image_data']
# Normalize and plot the data
image_data_norm = ImageNormalize(data['image_data'],
interval=ZScaleInterval(),
stretch=LinearStretch())
ax.imshow(data['image_data'],
norm=image_data_norm,
transform=ax.get_transform(data['wcs']),
cmap='gray')
ax.set_title(data['date'], fontsize=10)
# Add yellow circle around KBO position
c = Circle((data['ra'], data['dec']),
0.002,
edgecolor='yellow',
facecolor='none',
transform=ax.get_transform('world'))
ax.add_patch(c)
# Add grid and remove axis labels
ax.coords.grid(color='white', ls='solid')
ra, dec = ax.coords['ra'], ax.coords['dec']
ra.set_ticklabel_visible(False)
dec.set_ticklabel_visible(False)
ra.set_axislabel('')
dec.set_axislabel('')
except ValueError as ex:
print('Value Error: %s' % ex)
continue
axes = kbo_fig.get_axes()
axes[0].get_shared_x_axes().join(*axes)
axes[0].get_shared_y_axes().join(*axes)
plt.subplots_adjust(left=0.05, right=0.95, top=0.9)
plt.tight_layout(pad=0.5)
6. Putting it All Together¶
Now we can build the application. We can put the widgets together in a nice fashion using the HBox and VBox widgets then display everything.
from astropy.table import Table
range_slider = make_date_range_slider(cadc_table)
data_list = get_data_list(Table.from_pandas(
range_slider.children[0].selection),
ext=-1)
plt.ioff()
kbo_fig = kbo_figure(data_list, kbo_name)
plt.ion()
replot_button = make_replot_button(range_slider, kbo_fig, kbo_name)
# Organize the widgets
box_layout = widgets.Layout(align_items='stretch',
align_content='center',
border='none',
justify_content='space-around',
width='100%')
aperture_widget = make_aperture_text_widget()
silder_and_button = widgets.HBox([range_slider, replot_button],
layout=box_layout)
big_box = widgets.VBox([silder_and_button, kbo_fig.canvas, aperture_widget],
layout=box_layout)
display(big_box)
Use the interactive toolbar at the bottom of the plot to pan and zoom around. Note that each view is linked, so moving one axes view will move all the others. To select an area to preform an aperture sum, be sure to deselect the zoom/pan button then click in the centre of the area you want to analyze and drag out to the edge. The aperture sum label at the bottom should automatically be filled with the right value. Hopefully this tutorial has given you an understanding of how to build a simple astronomy application in a notebook.