Visualizing spectral cubes with yt

Extracting yt objects

The SpectralCube class includes a to_yt() method that makes is easy to return an object that can be used by yt to make volume renderings or other visualizations of the data. One common issue with volume rendering of spectral cubes is that you may not want pixels along the spectral axis to be given the same ‘3-d’ size as positional pixels, so the to_yt() method includes a spectral_factor argument that can be used to compress or expand the spectral axis.

The to_yt() method is used as follows:

>>> ytcube = cube.to_yt(spectral_factor=0.5)  
>>> ds = ytcube.dataset  

Warning

The API change in https://github.com/radio-astro-tools/spectral-cube/pull/129 affects the interpretation of the 0-pixel. There may be a 1-pixel offset between the yt cube and the SpectralCube

The ds object is then a yt object that can be used for rendering! By default the dataset is defined in pixel coordinates, going from 0.5 to n+0.5, as would be the case in ds9, for example. Along the spectral axis, this range will be modified if spectral_factor does not equal unity.

When working with datasets in yt, it may be useful to convert world coordinates to pixel coordinates, so that whenever you may have to input a position in yt (e.g., for slicing or volume rendering) you can get the pixel coordinate that corresponds to the desired world coordinate. For this purpose, the method world2yt() is provided:

>>> import astropy.units as u
>>> pix_coord = ytcube.world2yt([51.424522,
...                              30.723611,
...                              5205.18071],  # units of deg, deg, m/s
...                             )  

There is also a reverse method provided, yt2world():

>>> world_coord = ytcube.yt2world([ds.domain_center])  

which in this case would return the world coordinates of the center of the dataset in yt.

Note

The to_yt() method and its associated coordinate methods are compatible with both yt v. 2.x and v. 3.0 and following, but use of version 3.0 or later is recommended due to substantial improvements in support for FITS data. For more information on how yt handles FITS datasets, see the yt docs.

Visualization example

This section shows an example of a rendering script that can be used to produce a 3-d isocontour visualization using an object returned by to_yt():

import numpy as np
from spectral_cube import SpectralCube
import yt
import astropy.units as u

# Read in spectral cube
cube = SpectralCube.read('L1448_13CO.fits', format='fits')

# Extract the yt object from the SpectralCube instance
ytcube = cube.to_yt(spectral_factor=0.75)
ds = ytcube.dataset

# Set the number of levels, the minimum and maximum level and the width
# of the isocontours
n_v = 10
vmin = 0.05
vmax = 4.0
dv = 0.02

# Set up color transfer function
transfer = yt.ColorTransferFunction((vmin, vmax))
transfer.add_layers(n_v, dv, colormap='RdBu_r')

# Set up the camera parameters

# Derive the pixel coordinate of the desired center
# from the corresponding world coordinate
center = ytcube.world2yt([51.424522,
                          30.723611,
                          5205.18071])
direction = np.array([1.0, 0.0, 0.0])
width = 100.  # pixels
size = 1024

camera = ds.camera(center, direction, width, size, transfer,
                   fields=['flux'])

# Take a snapshot and save to a file
snapshot = camera.snapshot()
yt.write_bitmap(snapshot, 'cube_rendering.png', transpose=True)

You can move the camera around; see the yt camera docs.

Movie Making

There is a simple utility for quick movie making. The default movie is a rotation of the cube around one of the spatial axes, going from PP -> PV space and back.:

>>> cube = read('cube.fits', format='fits')  
>>> ytcube = cube.to_yt()  
>>> images = ytcube.quick_render_movie('outdir')  

The movie only does rotation, but it is a useful stepping-stone if you wish to learn how to use yt’s rendering system.

Example:

SketchFab Isosurface Contours

For data exploration, making movies can be tedious - it is difficult to control the camera and expensive to generate new renderings. Instead, creating a ‘model’ from the data and exporting that to SketchFab can be very useful. Only grayscale figures will be created with the quicklook code.

You need an account on sketchfab.com for this to work.:

>>> ytcube.quick_isocontour(title='GRS l=49 13CO 1 K contours', level=1.0)  

Here’s an example:

GRS l=49 13CO 1 K contours by keflavich on Sketchfab

You can also export locally to .ply and .obj files, which can be read by many programs (sketchfab, meshlab, blender). See the yt page for details.:

>>> ytcube.quick_isocontour(export_to='ply', filename='meshes.ply', level=1.0)  
>>> ytcube.quick_isocontour(export_to='obj', filename='meshes', level=1.0)