AAS234 WWT Workshop Notebook (June 2019; St. Loui, MO)

Please note: these notebooks are preserved for posterity but not kept up-to-date, so you might run into issues with older notebooks.

Data layers in pywwt.

In addition to the two and three-dimensional layer data included out of the box, pywwt also allows users to upload their own astropy Table objects and FITS image data for quick and seamless visualizations.

In this notebook, we'll use astroquery to download data in both formats and show how to plot them in multiple dimensions and against several all-sky, large-scale astronomical surveys.

Import and initialize the viewer.

(Select a cell and press ctrl + Enter to run it.)

In [ ]:
import time

from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.io import fits
from astropy.table import Table, Column
from astropy.time import Time, TimeDelta # note capitalization
from astropy.wcs import WCS

import numpy as np

from astroquery.nasa_exoplanet_archive import NasaExoplanetArchive
from astroquery.skyview import SkyView

from pywwt.jupyter import WWTJupyterWidget
In [ ]:
wwt = WWTJupyterWidget()
wwt

If you're using JupyterLab and not just a plain Jupyter notebook, you can move the WWT view to a separate window pane. This is extremely useful since it lets you keep on typing code without scrolling WWT out of view. Here's how you do that:

Right click and select "Create New View for Output"

If you don't get a menu when you right-click, or the menu doesn't look like the one pictured, you are using a plain Jupyter notebook and will have to scroll back and forth.

Plot K2’s footprint and exoplanet yield.

Here, we use pywwt to plot data from an astropy Table in both two and three dimensions.

To begin, we use a class from the astroquery package to query NASA's Exoplanet Archive and return an astropy Table that contains information about every confirmed exoplanet so far. It's a large table, so we'll only keep the columns necessary to run this example.

(If the astroquery call doesn't load, uncomment and run the cell following it to download the table from a pywwt-affiliated GitHub repository.)

In [ ]:
#xarch = NasaExoplanetArchive.get_confirmed_planets_table()
#xarch.keep_columns(['pl_hostname', 'pl_letter', 'st_dist', 'ra', 'dec'])
In [ ]:
xarch = Table.read('https://raw.githubusercontent.com/WorldWideTelescope/'
                   'pywwt-notebooks/master/aas-tutorials/data/k2_table.ecsv',
                   format='ascii.ecsv')

(Open a new cell above this one and see what the table looks like by entering xarch and running it.)

Next, we slice the table twice so it only contains entries of 'b' planets observed by K2.

K2 planets and stars are named using a consistent system:

K2 + 'system number' + 'letter'.

The letter for a star is typically 'a', the first exoplanet confirmed in the system is 'b', and so on.

We are plotting stars, but since the table only contains exoplanets (letters from 'b' onward), we will slice the table so only planets with the b suffix remain. That gives us one entry per system.

k2_pl searches the xarch column that lists star names and only keeps those with K2 at the beginning. Then, k2_st searches the suffix of the planet name and only keeps b planets.

(If you aren't familiar with list comprehensions, see Note 1 at the end of the notebook.)

In [ ]:
# select only k2 planets
k2_pl = xarch[ [row.startswith('K2') for row in xarch['pl_hostname']] ]

# select only 'b' planets from k2
k2_st = k2_pl[ ['b' in row for row in k2_pl['pl_letter']] ]

After that, we load the footprint of the K2 campaign fields on the sky using pywwt's add_fov method.

In [ ]:
field = wwt.add_fov(wwt.instruments.k2)

Allow a few seconds for the footprint to load, then watch as pywwtdraws it on the sky.

Once the footprint has loaded, we finally create the data layer that will plot our points in the viewer.

pywwt reads the locations from the keyword arguments lon_att and lat_att, and we feed it the names of the columns in k2_st that contain the stellar right ascensions and declinations. We can also set other attributes like color and size in the function call.

In [ ]:
lay = wwt.layers.add_table_layer(table=k2_st, frame='Sky',
                                lon_att='ra', lat_att='dec',
                                color='#c4d600', size_scale=40)

Do you see the points? If not, try making them larger.

In [ ]:
lay.size_scale *= 2

If you want to see the points and footprint without a sky background, pywwt also has a black layer available.

In [ ]:
wwt.foreground_opacity = 0
wwt.background = wwt.imagery.other.black

It is also possible to visualize data layers in pywwt's 3D solar system mode. To do so, we'll switch the view.

We will also remove the K2 footprint, as add_fov is designed to work primarily in 2D.

In [ ]:
field.remove()
wwt.set_view('solar system')

Zoom out, and you'll see K2's exoplanet-hosting stars in 3D. All of them are currently the same distance from Earth -- let's fix that.

Data layers also have attributes related to distance and altitude. We set alt_type to the correct type for our scenario, and let pywwt know which column of the table contains distance with alt_att.

In [ ]:
lay.alt_type = 'distance'
lay.alt_att = 'st_dist'

Zoom out again and enjoy the view of K2's stars scattered across the Milky Way at their proper distances from Earth.

Layer FITS images over existing all-sky surveys.

Here, we use pywwt to view a FITS image against the backdrop of multiple all-sky surveys.

We'll begin by resetting the viewer to its originial state. Then, since we're working with an infrared image, we'll switch the background to one of the same wavelength and make the foreground layer invisible for simplicity's sake.

In [ ]:
wwt.reset()
wwt.background = wwt.imagery.ir.twomass
wwt.foreground_opacity = 0

We again use astroquery to get our data, this time harnessing the SkyView service to retrieve a FITS cutout of a portion of M101 that experienced a Type Ia supernova in 2011.

The M101 cutout comes from a pre-event image from the highest wavelength band of 2MASS, an infrared all-sky survey. SkyView can also access of ther bands of 2MASS and many other surveys, as shown at the end of this document in Note 2.

(Again, if the astroquery call doesn't load, uncomment and run the cell following it to download the table from a pywwt-affiliated GitHub repository.)

In [ ]:
# Choose the size of the image in pixels
size = 500 # remember that larger sizes mean longer waits!

img_list = SkyView.get_images(position='SN 2011FE',
                              survey='2MASS-K', pixels=size)
In [ ]:
#img_list = [fits.open('https://github.com/WorldWideTelescope/'
#                      'pywwt-notebooks/blob/live-image-beta/'
#                      'tutorials/data/sn_2011fe_500.fits?raw=true')]

Because SkyView.get_images() returns a list, we have to index img_list once the download has finished in order to get to the FITS image data itself. (In this example, img_results is a list of just one.)

In [ ]:
sn11 = img_list[0]
sn11.info()

In the most basic cases, this is all that's needed to create an image layer.

If you have the PrimaryHDU, or "Primary Header Data Unit", pywwt only needs one command to plot it. Once uploaded, the viewer will take you to the proper location and scale to see your image.

In [ ]:
img = wwt.layers.add_image_layer(sn11)

(It's also possible to upload FITS images from a local file by using image='the/path/to/your/file' as the method's argument instead of an existing Python object name. pywwt will then automatically do the examination of the header and image data that we did in previous cells.)

We can also separate PrimaryHDU into its header and image data components and make edits to either before creating a layer. We'll explore that by splitting them into different variables and deleting our first image layer from the viewer.

In [ ]:
sn11_data = sn11[0].data.copy()
sn11_header = sn11[0].header.copy()
img.remove()

pywwt requires the image data to be structured as a 2-D array, which is what we see below. The dimensions of our array are (size x size), just as we specified earlier.

In [ ]:
print(sn11_data.shape)
sn11_data

We can perform operations on this array like scrubbing a column of pixels to, say, censor out an artifact that we don't want to display.

In [ ]:
sn11_data[0:size, 50:100] = np.median(sn11_data)

The header contains base-level data about the image such as its coordinates and the projection system it follows. We will need it to create an astropy.wcs.WCS object, which pywwt uses to convert the image to its preferred coordinate system and orient the data properly on its 2-D sky projection.

In [ ]:
sn11_header

Speaking of which, let's create the WCS (World Coordinate System) object now. As stated before, all we need is the header data; astropy's internals do the rest of the heavy lifting.

In [ ]:
WCS(sn11_header)

We can then bring these separate header and data objects together in a tuple to run the add_image_layer method again. You should notice a vertical band that wasn't there before -- we created that by manipulating the image data array.

In [ ]:
img = wwt.layers.add_image_layer(image=(sn11_data, WCS(sn11_header)))

Once the image has loaded, we are able to edit a few of its attributes: opacity, vmin/vmax (the minimum and maximum flux values covered by the colormap), and stretch (the scaling function used to display the data). There are more options coming as we continue development.

Changing opacity shows us that our image is indeed aligned with the background survey.

In [ ]:
img.opacity = .4
In [ ]:
img.opacity = .6

Shifting vmin and vmax changes the lower and upper bounds of our grayscale colormap.

pywwt automatically sets these attributes such that ~99% of points in the image data array are covered, but we can adjust them depending on what exactly we'd like to emphasize.

In [ ]:
print(np.percentile(sn11_data, .5), np.percentile(sn11_data, 99.5))
print(img.vmin, img.vmax)
In [ ]:
img.vmin = 360
In [ ]:
img.vmax = 370

Valid stretch values are: linear, log, power, sqrt, and histeq. pywwt chooses a linear stretch by default, but again, we can adjust it to fit our individual needs -- aesthetic, scientific, or anywhere in between.

In [ ]:
img.stretch = 'log'

The great advantage of being able to upload personal FITS images into pywwt is that we can easily compare them with multiple surveys at once.

Since our FITS image is in infrared, we can use a background and a foreground all-sky survey from pywwt to see M101 in three wavelengths at once. (Please give them a few seconds to load.)

In [ ]:
wwt.background = wwt.imagery.visible.sdss
wwt.foreground = wwt.imagery.gamma.fermi

To get the right balance between the all-sky layers, we can manually change the opacity of the foreground...

In [ ]:
wwt.foreground_opacity = .7

...or we can use a Jupyter-exclusive option to do the same -- and more -- using live sliders and menus.

In [ ]:
wwt.layer_controls

We can use pywwt's ability to add annotations to the viewer and SkyCoord's location functionality to visually indicate the spot where our supernova took place.

In [ ]:
sn11_coord = SkyCoord.from_name('Sn 2011FE')
circle = wwt.add_circle(sn11_coord, radius=.003*u.deg)

As in the K2 example, we're able to change attributes of the annotation after its creation and remove it when we no longer need it.

In [ ]:
circle.line_width = 6 * u.pixel
circle.line_color = '#e56020'
circle.radius = wwt.get_fov() / 20
circle.set_center(wwt.get_center())
In [ ]:
circle.remove()

Finally, remember that pywwt is four-dimensional visualization engine -- this means we can also visualize the progression of time.

Using astropy Time objects to change the radius of our circle over time, we can make a basic animation of how SN 2011fe looked in the months immediately surrounding its detonation and view it in a matter of seconds.

(This animation is qualitative and is only meant to approximate how the supernova looked.)

In [ ]:
circ = wwt.add_circle(sn11_coord, fill=True, radius=.003*u.deg)
In [ ]:
# Save benchmark times for animation
start = Time('2011-08-01', format='iso')
peak = Time('2011-09-13', format='iso')
finish = Time('2011-10-26', format='iso')
span = peak - start

now = start
step = (finish - start) / 20  - TimeDelta(.01 * u.second)
In [ ]:
# Loop through a few months in a few seconds of real time
max_rad = circ.radius

while now <= finish:
    wwt.set_current_time(now)
    
    circ.radius = max_rad * (1 - np.abs(now - peak) / span)
    time.sleep(.1)
    
    now += step
circ.remove()

Is it working? If so, congratulations, enjoy the view, and stay in touch! Find us at the AAS booth, sign up for our mailing list, discuss with other users at the WWT forum, and contribute suggestions and fixes at our GitHub repository.

Note 1

The constructions used to build k2_pl and k2_st are called "list comprehensions," and they are shortcuts for building lists. The regular, equivalent method to create k2_st, looks like this:

index_list = []
for i, row in enumerate(xarch, 0):
    if row['pl_hostname'].startswith('K2'):
        index_list.append(i)

k2_pl = xarch[index_list]

Note 2

You can see all of the layers available layers in SkyView by running SkyView.list_surveys() or SkyView.survey_dict(), but be aware that these will take some time to load.

To get links to image files from all bands of a particular survey, run:

SkyView.get_image_list('Sn 2011FE', survey=SkyView.survey_dict['IR:2MASS class='])

You can substitute the dictionary entry for whichever surveys are present in SkyView.survey_dict(). For example, if you wanted images from all bands of SDSS, the survey keyword argument in the call above would change to SkyView.survey_dict['Optical:SDSS']

Credits

This notebook was prepared by O. Justin Otor.