A Few Complex Plots

The built-in plotting functionality covers the very simple use cases that are most common. These scripts will demonstrate how to construct more complex plots or publication-quality plots. In many cases these show how to make multi-panel plots.

Multi-Width Image

This is a simple recipe to show how to open a dataset and then plot slices through it at varying widths. See Slice Plots for more information.

(multi_width_image.py)

import yt

# Load the dataset.
ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")

# Create a slice plot for the dataset.  With no additional arguments,
# the width will be the size of the domain and the center will be the
# center of the simulation box
slc = yt.SlicePlot(ds, 'z', 'density')

# Create a list of a couple of widths and units.
# (N.B. Mpc (megaparsec) != mpc (milliparsec)
widths = [(1, 'Mpc'),
          (15, 'kpc')]

# Loop through the list of widths and units.
for width, unit in widths:

    # Set the width.
    slc.set_width(width, unit)

    # Write out the image with a unique name.
    slc.save("%s_%010d_%s" % (ds, width, unit))

zoomFactors = [2,4,5]

# recreate the original slice
slc = yt.SlicePlot(ds, 'z', 'density')

for zoomFactor in zoomFactors:

    # zoom in
    slc.zoom(zoomFactor)

    # Write out the image with a unique name.
    slc.save("%s_%i" % (ds, zoomFactor))
../_images/multi_width_image__galaxy0030_0000000001_Mpc_Slice_z_density.png ../_images/multi_width_image__galaxy0030_0000000015_kpc_Slice_z_density.png ../_images/multi_width_image__galaxy0030_2_Slice_z_density.png ../_images/multi_width_image__galaxy0030_4_Slice_z_density.png ../_images/multi_width_image__galaxy0030_5_Slice_z_density.png

Multipanel with Axes Labels

This illustrates how to use a SlicePlot to control a multipanel plot. This plot uses axes labels to illustrate the length scales in the plot. See Slice Plots and the Matplotlib AxesGrid Object for more information.

(multiplot_2x2.py)

import yt
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import AxesGrid

fn = "IsolatedGalaxy/galaxy0030/galaxy0030"
ds = yt.load(fn) # load data

fig = plt.figure()

# See http://matplotlib.org/mpl_toolkits/axes_grid/api/axes_grid_api.html
# These choices of keyword arguments produce a four panel plot that includes
# four narrow colorbars, one for each plot.  Axes labels are only drawn on the
# bottom left hand plot to avoid repeating information and make the plot less
# cluttered.
grid = AxesGrid(fig, (0.075,0.075,0.85,0.85),
                nrows_ncols = (2, 2),
                axes_pad = 1.0,
                label_mode = "1",
                share_all = True,
                cbar_location="right",
                cbar_mode="each",
                cbar_size="3%",
                cbar_pad="0%")

fields = ['density', 'velocity_x', 'velocity_y', 'velocity_magnitude']

# Create the plot.  Since SlicePlot accepts a list of fields, we need only
# do this once.
p = yt.SlicePlot(ds, 'z', fields)

# Velocity is going to be both positive and negative, so let's make these
# slices use a linear colorbar scale
p.set_log('velocity_x', False)
p.set_log('velocity_y', False)

p.zoom(2)

# For each plotted field, force the SlicePlot to redraw itself onto the AxesGrid
# axes.
for i, field in enumerate(fields):
    plot = p.plots[field]
    plot.figure = fig
    plot.axes = grid[i].axes
    plot.cax = grid.cbar_axes[i]

# Finally, redraw the plot on the AxesGrid axes.
p._setup_plots()

plt.savefig('multiplot_2x2.png')
../_images/multiplot_2x2__multiplot_2x2.png

Time Series Multipanel

This illustrates how to create a multipanel plot of a time series dataset. See Projection Plots, Time Series Analysis, and the Matplotlib AxesGrid Object for more information.

(multiplot_2x2_time_series.py)

import yt
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import AxesGrid

fns = ['Enzo_64/DD0005/data0005','Enzo_64/DD0015/data0015', 'Enzo_64/DD0025/data0025', 'Enzo_64/DD0035/data0035']

fig = plt.figure()

# See http://matplotlib.org/mpl_toolkits/axes_grid/api/axes_grid_api.html
# These choices of keyword arguments produce a four panel plot with a single
# shared narrow colorbar on the right hand side of the multipanel plot. Axes
# labels are drawn for all plots since we're slicing along different directions
# for each plot.
grid = AxesGrid(fig, (0.075,0.075,0.85,0.85),
                nrows_ncols = (2, 2),
                axes_pad = 0.05,
                label_mode = "L",
                share_all = True,
                cbar_location="right",
                cbar_mode="single",
                cbar_size="3%",
                cbar_pad="0%")

for i, fn in enumerate(fns):
    # Load the data and create a single plot
    ds = yt.load(fn) # load data
    p = yt.ProjectionPlot(ds, 'z', 'density', width=(55, 'Mpccm'))

    # Ensure the colorbar limits match for all plots
    p.set_zlim('density', 1e-4, 1e-2)

    # This forces the ProjectionPlot to redraw itself on the AxesGrid axes.
    plot = p.plots['density']
    plot.figure = fig
    plot.axes = grid[i].axes
    plot.cax = grid.cbar_axes[i]

    # Finally, this actually redraws the plot.
    p._setup_plots()

plt.savefig('multiplot_2x2_time_series.png')
../_images/multiplot_2x2_time_series__multiplot_2x2_time_series.png

Mutiple Slice Multipanel

This illustrates how to create a multipanel plot of slices along the coordinate axes. To focus on what’s happening in the x-y plane, we make an additional Temperature slice for the bottom-right subpanel. See Slice Plots and the Matplotlib AxesGrid Object for more information.

(multiplot_2x2_coordaxes_slice.py)

import yt
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import AxesGrid

fn = "IsolatedGalaxy/galaxy0030/galaxy0030"
ds = yt.load(fn) # load data

fig = plt.figure()

# See http://matplotlib.org/mpl_toolkits/axes_grid/api/axes_grid_api.html
# These choices of keyword arguments produce two colorbars, both drawn on the
# right hand side.  This means there are only two colorbar axes, one for Density
# and another for temperature.  In addition, axes labels will be drawn for all
# plots.
grid = AxesGrid(fig, (0.075,0.075,0.85,0.85),
                nrows_ncols = (2, 2),
                axes_pad = 1.0,
                label_mode = "all",
                share_all = True,
                cbar_location="right",
                cbar_mode="edge",
                cbar_size="5%",
                cbar_pad="0%")

cuts = ['x', 'y', 'z', 'z']
fields = ['density', 'density', 'density', 'temperature']

for i, (direction, field) in enumerate(zip(cuts, fields)):
    # Load the data and create a single plot
    p = yt.SlicePlot(ds, direction, field)
    p.zoom(40)

    # This forces the ProjectionPlot to redraw itself on the AxesGrid axes.
    plot = p.plots[field]
    plot.figure = fig
    plot.axes = grid[i].axes

    # Since there are only two colorbar axes, we need to make sure we don't try
    # to set the temperature colorbar to cbar_axes[4], which would if we used i
    # to index cbar_axes, yielding a plot without a temperature colorbar.
    # This unecessarily redraws the Density colorbar three times, but that has
    # no effect on the final plot.
    if field == 'density':
        plot.cax = grid.cbar_axes[0]
    elif field == 'temperature':
        plot.cax = grid.cbar_axes[1]

    # Finally, redraw the plot.
    p._setup_plots()

plt.savefig('multiplot_2x2_coordaxes_slice.png')
../_images/multiplot_2x2_coordaxes_slice__multiplot_2x2_coordaxes_slice.png

Multi-Plot Slice and Projections

This shows how to combine multiple slices and projections into a single image, with detailed control over colorbars, titles and color limits. See Slice Plots and Projection Plots for more information.

(multi_plot_slice_and_proj.py)

import yt
import numpy as np
from yt.visualization.base_plot_types import get_multi_plot
import matplotlib.colorbar as cb
from matplotlib.colors import LogNorm

fn = "GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150" # dataset to load
orient = 'horizontal'

ds = yt.load(fn) # load data

# There's a lot in here:
#   From this we get a containing figure, a list-of-lists of axes into which we
#   can place plots, and some axes that we'll put colorbars.
# We feed it:
#   Number of plots on the x-axis, number of plots on the y-axis, and how we
#   want our colorbars oriented.  (This governs where they will go, too.
#   bw is the base-width in inches, but 4 is about right for most cases.
fig, axes, colorbars = get_multi_plot(3, 2, colorbar=orient, bw = 4)

slc = yt.SlicePlot(ds, 'z', fields=["density","temperature","velocity_magnitude"])
proj = yt.ProjectionPlot(ds, 'z', "density", weight_field="density")

slc_frb = slc.data_source.to_frb((1.0, "Mpc"), 512)
proj_frb = proj.data_source.to_frb((1.0, "Mpc"), 512)

dens_axes = [axes[0][0], axes[1][0]]
temp_axes = [axes[0][1], axes[1][1]]
vels_axes = [axes[0][2], axes[1][2]]

for dax, tax, vax in zip(dens_axes, temp_axes, vels_axes) :

    dax.xaxis.set_visible(False)
    dax.yaxis.set_visible(False)
    tax.xaxis.set_visible(False)
    tax.yaxis.set_visible(False)
    vax.xaxis.set_visible(False)
    vax.yaxis.set_visible(False)

# Converting our Fixed Resolution Buffers to numpy arrays so that matplotlib
# can render them

slc_dens = np.array(slc_frb['density'])
proj_dens = np.array(proj_frb['density'])
slc_temp = np.array(slc_frb['temperature'])
proj_temp = np.array(proj_frb['temperature'])
slc_vel = np.array(slc_frb['velocity_magnitude'])
proj_vel = np.array(proj_frb['velocity_magnitude'])

plots = [dens_axes[0].imshow(slc_dens, origin='lower', norm=LogNorm()),
         dens_axes[1].imshow(proj_dens, origin='lower', norm=LogNorm()),
         temp_axes[0].imshow(slc_temp, origin='lower'),
         temp_axes[1].imshow(proj_temp, origin='lower'),
         vels_axes[0].imshow(slc_vel, origin='lower', norm=LogNorm()),
         vels_axes[1].imshow(proj_vel, origin='lower', norm=LogNorm())]

plots[0].set_clim((1.0e-27,1.0e-25))
plots[0].set_cmap("bds_highcontrast")
plots[1].set_clim((1.0e-27,1.0e-25))
plots[1].set_cmap("bds_highcontrast")
plots[2].set_clim((1.0e7,1.0e8))
plots[2].set_cmap("hot")
plots[3].set_clim((1.0e7,1.0e8))
plots[3].set_cmap("hot")
plots[4].set_clim((1e6, 1e8))
plots[4].set_cmap("gist_rainbow")
plots[5].set_clim((1e6, 1e8))
plots[5].set_cmap("gist_rainbow")

titles=[r'$\mathrm{Density}\ (\mathrm{g\ cm^{-3}})$',
        r'$\mathrm{Temperature}\ (\mathrm{K})$',
        r'$\mathrm{Velocity Magnitude}\ (\mathrm{cm\ s^{-1}})$']

for p, cax, t in zip(plots[0:6:2], colorbars, titles):
    cbar = fig.colorbar(p, cax=cax, orientation=orient)
    cbar.set_label(t)

# And now we're done!
fig.savefig("%s_3x2" % ds)
../_images/multi_plot_slice_and_proj__sloshing_nomag2_hdf5_plt_cnt_0150_3x2.png

Advanced Multi-Plot Multi-Panel

This produces a series of slices of multiple fields with different color maps and zlimits, and makes use of the FixedResolutionBuffer. While this is more complex than the equivalent plot collection-based solution, it allows for a lot more flexibility. Every part of the script uses matplotlib commands, allowing its full power to be exercised. See Slice Plots and Projection Plots for more information.

(multi_plot_3x2_FRB.py)

import yt
import numpy as np
from yt.visualization.api import get_multi_plot
import matplotlib.colorbar as cb
from matplotlib.colors import LogNorm

fn = "Enzo_64/RD0006/RedshiftOutput0006" # dataset to load

# load data and get center value and center location as maximum density location
ds = yt.load(fn)
v, c = ds.find_max("density")

# set up our Fixed Resolution Buffer parameters: a width, resolution, and center
width = (1.0, 'unitary')
res = [1000, 1000]
#  get_multi_plot returns a containing figure, a list-of-lists of axes
#   into which we can place plots, and some axes that we'll put
#   colorbars.

#  it accepts: # of x-axis plots, # of y-axis plots, and how the
#  colorbars are oriented (this also determines where they go: below
#  in the case of 'horizontal', on the right in the case of
#  'vertical'), bw is the base-width in inches (4 is about right for
#  most cases)

orient = 'horizontal'
fig, axes, colorbars = get_multi_plot( 2, 3, colorbar=orient, bw = 6)

# Now we follow the method of "multi_plot.py" but we're going to iterate
# over the columns, which will become axes of slicing.
plots = []
for ax in range(3):
    sli = ds.slice(ax, c[ax])
    frb = sli.to_frb(width, res)
    den_axis = axes[ax][0]
    temp_axis = axes[ax][1]

    # here, we turn off the axes labels and ticks, but you could
    # customize further.
    for ax in (den_axis, temp_axis):
        ax.xaxis.set_visible(False)
        ax.yaxis.set_visible(False)

    # converting our fixed resolution buffers to NDarray so matplotlib can
    # render them
    dens = np.array(frb['density'])
    temp = np.array(frb['temperature'])

    plots.append(den_axis.imshow(dens, norm=LogNorm()))
    plots[-1].set_clim((5e-32, 1e-29))
    plots[-1].set_cmap("bds_highcontrast")

    plots.append(temp_axis.imshow(temp, norm=LogNorm()))
    plots[-1].set_clim((1e3, 1e8))
    plots[-1].set_cmap("hot")

# Each 'cax' is a colorbar-container, into which we'll put a colorbar.
# the zip command creates triples from each element of the three lists
# .  Note that it cuts off after the shortest iterator is exhausted,
# in this case, titles.
titles=[r'$\mathrm{Density}\ (\mathrm{g\ cm^{-3}})$', r'$\mathrm{temperature}\ (\mathrm{K})$']
for p, cax, t in zip(plots, colorbars,titles):
    # Now we make a colorbar, using the 'image' we stored in plots
    # above. note this is what is *returned* by the imshow method of
    # the plots.
    cbar = fig.colorbar(p, cax=cax, orientation=orient)
    cbar.set_label(t)

# And now we're done!
fig.savefig("%s_3x2.png" % ds)
../_images/multi_plot_3x2_FRB__RedshiftOutput0006_3x2.png

Off-Axis Projection (an alternate method)

This recipe demonstrates how to take an image-plane line integral along an arbitrary axis in a simulation. This uses alternate machinery than the standard PlotWindow interface to create an off-axis projection as demonstrated in this recipe.

(offaxis_projection.py)

import yt
import numpy as np

# Load the dataset.
ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")

# Choose a center for the render.
c = [0.5, 0.5, 0.5]

# Our image plane will be normal to some vector.  For things like collapsing
# objects, you could set it the way you would a cutting plane -- but for this
# dataset, we'll just choose an off-axis value at random.  This gets normalized
# automatically.
L = [1.0, 0.0, 0.0]

# Our "width" is the width of the image plane as well as the depth.
# The first element is the left to right width, the second is the
# top-bottom width, and the last element is the back-to-front width
# (all in code units)
W = [0.04,0.04,0.4]

# The number of pixels along one side of the image.
# The final image will have Npixel^2 pixels.
Npixels = 512

# Create the off axis projection.
# Setting no_ghost to False speeds up the process, but makes a
# slighly lower quality image.
image = yt.off_axis_projection(ds, c, L, W, Npixels, "density", no_ghost=False)

# Write out the final image and give it a name
# relating to what our dataset is called.
# We save the log of the values so that the colors do not span
# many orders of magnitude.  Try it without and see what happens.
yt.write_image(np.log10(image), "%s_offaxis_projection.png" % ds)
../_images/offaxis_projection__galaxy0030_offaxis_projection.png

Off-Axis Projection with a Colorbar (an alternate method)

This recipe shows how to generate a colorbar with a projection of a dataset from an arbitrary projection angle (so you are not confined to the x, y, and z axes).

This uses alternate machinery than the standard PlotWindow interface to create an off-axis projection as demonstrated in this recipe.

(offaxis_projection_colorbar.py)

import yt
import numpy as np

fn = "IsolatedGalaxy/galaxy0030/galaxy0030" # dataset to load

ds = yt.load(fn) # load data

# Now we need a center of our volume to render.  Here we'll just use
# 0.5,0.5,0.5, because volume renderings are not periodic.
c = [0.5, 0.5, 0.5]

# Our image plane will be normal to some vector.  For things like collapsing
# objects, you could set it the way you would a cutting plane -- but for this
# dataset, we'll just choose an off-axis value at random.  This gets normalized
# automatically.
L = [0.5, 0.4, 0.7]

# Our "width" is the width of the image plane as well as the depth.
# The first element is the left to right width, the second is the
# top-bottom width, and the last element is the back-to-front width
# (all in code units)
W = [0.04,0.04,0.4]

# The number of pixels along one side of the image.
# The final image will have Npixel^2 pixels.
Npixels = 512

# Now we call the off_axis_projection function, which handles the rest.
# Note that we set no_ghost equal to False, so that we *do* include ghost
# zones in our data.  This takes longer to calculate, but the results look
# much cleaner than when you ignore the ghost zones.
# Also note that we set the field which we want to project as "density", but
# really we could use any arbitrary field like "temperature", "metallicity"
# or whatever.
image = yt.off_axis_projection(ds, c, L, W, Npixels, "density", no_ghost=False)

# Image is now an NxN array representing the intensities of the various pixels.
# And now, we call our direct image saver.  We save the log of the result.
yt.write_projection(image, "offaxis_projection_colorbar.png",
                    colorbar_label="Column Density (cm$^{-2}$)")
../_images/offaxis_projection_colorbar__offaxis_projection_colorbar.png

Thin-Slice Projections

This recipe is an example of how to project through only a given data object, in this case a thin region, and then display the result. See Projection Plots and Available Objects for more information.

(thin_slice_projection.py)

import yt

# Load the dataset.
ds = yt.load("Enzo_64/DD0030/data0030")

# Make a projection that is the full width of the domain,
# but only 5 Mpc in depth.  This is done by creating a
# region object with this exact geometry and providing it
# as a data_source for the projection.

# Center on the domain center
center = ds.domain_center

# First make the left and right corner of the region based
# on the full domain.
left_corner = ds.domain_left_edge
right_corner = ds.domain_right_edge

# Now adjust the size of the region along the line of sight (x axis).
depth = ds.quan(5.0,'Mpc')
left_corner[0] = center[0] - 0.5 * depth
right_corner[0] = center[0] + 0.5 * depth

# Create the region
region = ds.box(left_corner, right_corner)

# Create a density projection and supply the region we have just created.
# Only cells within the region will be included in the projection.
# Try with another data container, like a sphere or disk.
plot = yt.ProjectionPlot(ds, "x", "density", weight_field="density",
                         data_source=region)

# Save the image with the keyword.
plot.save("Thin_Slice")
../_images/thin_slice_projection__Thin_Slice_Projection_x_density_density.png

Plotting Particles Over Fluids

This recipe demonstrates how to overplot particles on top of a fluid image. See Overplotting Particle Positions for more information.

(overplot_particles.py)

import yt

# Load the dataset.
ds = yt.load("Enzo_64/DD0043/data0043")

# Make a density projection centered on the 'm'aximum density location
# with a width of 10 Mpc..
p = yt.ProjectionPlot(ds, "y", "density", center='m', width=(10, 'Mpc'))

# Modify the projection
# The argument specifies the region along the line of sight
# for which particles will be gathered.
# 1.0 signifies the entire domain in the line of sight
# p.annotate_particles(1.0)
# but in this case we only go 10 Mpc in depth
p.annotate_particles((10, 'Mpc'))

# Save the image.
# Optionally, give a string as an argument
# to name files with a keyword.
p.save()
../_images/overplot_particles__data0043_Projection_y_density.png

Plotting Grid Edges Over Fluids

This recipe demonstrates how to overplot grid boxes on top of a fluid image. Each level is represented with a different color from white (low refinement) to black (high refinement). One can change the colormap used for the grids colors by using the cmap keyword (or set it to None to get all grid edges as black). See Overplot Grids for more information.

(overplot_grids.py)

import yt

# Load the dataset.
ds = yt.load("Enzo_64/DD0043/data0043")

# Make a density projection.
p = yt.ProjectionPlot(ds, "y", "density")

# Modify the projection
# The argument specifies the region along the line of sight
# for which particles will be gathered.
# 1.0 signifies the entire domain in the line of sight.
p.annotate_grids()

# Save the image.
# Optionally, give a string as an argument
# to name files with a keyword.
p.save()
../_images/overplot_grids__data0043_Projection_y_density.png

Overplotting Velocity Vectors

This recipe demonstrates how to plot velocity vectors on top of a slice. See Overplot Quivers for the Velocity Field for more information.

(velocity_vectors_on_slice.py)

import yt

# Load the dataset.
ds = yt.load("GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150")

p = yt.SlicePlot(ds, "x", "density")

# Draw a velocity vector every 16 pixels.
p.annotate_velocity(factor = 16)
p.save()
../_images/velocity_vectors_on_slice__sloshing_nomag2_hdf5_plt_cnt_0150_Slice_x_density.png

Overplotting Contours

This is a simple recipe to show how to open a dataset, plot a slice through it, and add contours of another quantity on top. See Overplot Contours for more information.

(contours_on_slice.py)

import yt

# first add density contours on a density slice
ds = yt.load("GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150")

# add density contours on the density slice.
p = yt.SlicePlot(ds, "x", "density")
p.annotate_contour("density")
p.save()

# then add temperature contours on the same density slice
p = yt.SlicePlot(ds, "x", "density")
p.annotate_contour("temperature")
p.save(str(ds)+'_T_contour')
../_images/contours_on_slice__sloshing_nomag2_hdf5_plt_cnt_0150_Slice_x_density.png ../_images/contours_on_slice__sloshing_nomag2_hdf5_plt_cnt_0150_T_contour_Slice_x_density.png

Simple Contours in a Slice

Sometimes it is useful to plot just a few contours of a quantity in a dataset. This shows how one does this by first making a slice, adding contours, and then hiding the colormap plot of the slice to leave the plot containing only the contours that one has added. See Overplot Contours for more information.

(simple_contour_in_slice.py)

import yt

# Load the data file.
ds = yt.load("Sedov_3d/sedov_hdf5_chk_0002")

# Make a traditional slice plot.
sp = yt.SlicePlot(ds, "x", "density")

# Overlay the slice plot with thick red contours of density.
sp.annotate_contour("density", ncont=3, clim=(1e-2,1e-1), label=True,
                    plot_args={"colors": "red",
                               "linewidths": 2})

# What about some nice temperature contours in blue?
sp.annotate_contour("temperature", ncont=3, clim=(1e-8,1e-6), label=True,
                    plot_args={"colors": "blue",
                               "linewidths": 2})

# This is the plot object.
po = sp.plots["density"]

# Turn off the colormap image, leaving just the contours.
po.axes.images[0].set_visible(False)

# Remove the colorbar and its label.
po.figure.delaxes(po.figure.axes[1])

# Save it and ask for a close fit to get rid of the space used by the colorbar.
sp.save(mpl_kwargs={'bbox_inches':'tight'})
../_images/simple_contour_in_slice__sedov_hdf5_chk_0002_Slice_x_density.png

Styling Radial Profile Plots

This recipe demonstrates a method of calculating radial profiles for several quantities, styling them and saving out the resultant plot. See 1D Profile Plots for more information.

(radial_profile_styles.py)

import yt
import matplotlib.pyplot as plt

ds = yt.load("GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150")

# Get a sphere object

sp = ds.sphere(ds.domain_center, (500., "kpc"))

# Bin up the data from the sphere into a radial profile

rp = yt.create_profile(sp, 'radius', ['density', 'temperature'],
                       units = {'radius': 'kpc'},
                       logs = {'radius': False})

# Make plots using matplotlib

fig = plt.figure()
ax = fig.add_subplot(111)

# Plot the density as a log-log plot using the default settings
dens_plot = ax.loglog(rp.x.value, rp["density"].value)

# Here we set the labels of the plot axes

ax.set_xlabel(r"$\mathrm{r\ (kpc)}$")
ax.set_ylabel(r"$\mathrm{\rho\ (g\ cm^{-3})}$")

# Save the default plot

fig.savefig("density_profile_default.png" % ds)

# The "dens_plot" object is a list of plot objects. In our case we only have one,
# so we index the list by '0' to get it.

# Plot using dashed red lines

dens_plot[0].set_linestyle("--")
dens_plot[0].set_color("red")

fig.savefig("density_profile_dashed_red.png")

# Increase the line width and add points in the shape of x's

dens_plot[0].set_linewidth(5)
dens_plot[0].set_marker("x")
dens_plot[0].set_markersize(10)

fig.savefig("density_profile_thick_with_xs.png")
../_images/radial_profile_styles__density_profile_dashed_red.png ../_images/radial_profile_styles__density_profile_default.png ../_images/radial_profile_styles__density_profile_thick_with_xs.png

Customized Profile Plot

This recipe demonstrates how to create a fully customized 1D profile object using the create_profile() function and then create a ProfilePlot using the customized profile. This illustrates how a ProfilePlot created this way inherits the properties of the profile it is constructed from. See 1D Profile Plots for more information.

(customized_profile_plot.py)

import yt
import yt.units as u

ds = yt.load('HiresIsolatedGalaxy/DD0044/DD0044')

center = [0.53, 0.53, 0.53]
normal = [0,0,1]
radius = 40*u.kpc
height = 5*u.kpc

disk = ds.disk(center, [0,0,1], radius, height)

profile = yt.create_profile(
    data_source=disk,
    bin_fields=["radius"],
    fields=["cylindrical_tangential_velocity_absolute"],
    n_bins=256,
    units=dict(radius="kpc",
               cylindrical_tangential_velocity_absolute="km/s"),
    logs=dict(radius=False),
    weight_field='cell_mass',
    extrema=dict(radius=(0,40)),
    )

plot = yt.ProfilePlot.from_profiles(profile)

plot.set_log('cylindrical_tangential_velocity_absolute', False)
plot.set_ylim('cylindrical_tangential_velocity_absolute', 60, 160)

plot.save()
../_images/customized_profile_plot__DD0044_1d-Profile_radius_cylindrical_tangential_velocity_absolute.png

Customized Phase Plot

Similar to the recipe above, this demonstrates how to create a fully customized 2D profile object using the create_profile() function and then create a PhasePlot using the customized profile object. This illustrates how a PhasePlot created this way inherits the properties of the profile object from which it is constructed. See 2D Phase Plots for more information.

(customized_phase_plot.py)

import yt
import yt.units as u

ds = yt.load('HiresIsolatedGalaxy/DD0044/DD0044')

center = [0.53, 0.53, 0.53]
normal = [0,0,1]
radius = 40*u.kpc
height = 2*u.kpc

disk = ds.disk(center, [0,0,1], radius, height)

profile = yt.create_profile(
    data_source=disk,
    bin_fields=["radius", "cylindrical_tangential_velocity"],
    fields=["cell_mass"],
    n_bins=256,
    units=dict(radius="kpc",
               cylindrical_tangential_velocity="km/s",
               cell_mass="Msun"),
    logs=dict(radius=False,
              cylindrical_tangential_velocity=False),
    weight_field=None,
    extrema=dict(radius=(0,40),
                 cylindrical_tangential_velocity=(-250, 250)),
    )

plot = yt.PhasePlot.from_profile(profile)
plot.set_cmap("cell_mass", "YlOrRd")

plot.save()
../_images/customized_phase_plot__DD0044_2d-Profile_radius_cylindrical_tangential_velocity_cell_mass.png

Moving a Volume Rendering Camera

In this recipe, we move a camera through a domain and take multiple volume rendering snapshots. This recipe uses an unstructured mesh dataset (see Unstructured Mesh Rendering), which makes it easier to visualize what the Camera is doing, but you can manipulate the Camera for other dataset types in exactly the same manner.

See Moving and Orienting the Camera for more information.

(camera_movement.py)

import yt
import numpy as np

ds = yt.load("MOOSE_sample_data/out.e-s010")
sc = yt.create_scene(ds)
cam = sc.camera

# save an image at the starting position
frame = 0
sc.save('camera_movement_%04i.png' % frame)
frame += 1

# Zoom out by a factor of 2 over 5 frames
for _ in cam.iter_zoom(0.5, 5):
    sc.render()
    sc.save('camera_movement_%04i.png' % frame)
    frame += 1

# Move to the position [-10.0, 10.0, -10.0] over 5 frames
pos = ds.arr([-10.0, 10.0, -10.0], 'code_length')
for _ in cam.iter_move(pos, 5):
    sc.render()
    sc.save('camera_movement_%04i.png' % frame)
    frame += 1

# Rotate by 180 degrees over 5 frames
for _ in cam.iter_rotate(np.pi, 5):
    sc.render()
    sc.save('camera_movement_%04i.png' % frame)
    frame += 1
../_images/camera_movement__camera_movement_0000.png ../_images/camera_movement__camera_movement_0001.png ../_images/camera_movement__camera_movement_0002.png ../_images/camera_movement__camera_movement_0003.png ../_images/camera_movement__camera_movement_0004.png ../_images/camera_movement__camera_movement_0005.png ../_images/camera_movement__camera_movement_0006.png ../_images/camera_movement__camera_movement_0007.png ../_images/camera_movement__camera_movement_0008.png ../_images/camera_movement__camera_movement_0009.png ../_images/camera_movement__camera_movement_0010.png ../_images/camera_movement__camera_movement_0011.png ../_images/camera_movement__camera_movement_0012.png ../_images/camera_movement__camera_movement_0013.png ../_images/camera_movement__camera_movement_0014.png ../_images/camera_movement__camera_movement_0015.png

Volume Rendering with Custom Camera

In this recipe we modify the Simple Volume Rendering recipe to use customized camera properties. See 3D Visualization and Volume Rendering for more information.

(custom_camera_volume_rendering.py)

import yt

# Load the dataset
ds = yt.load("Enzo_64/DD0043/data0043")

# Create a volume rendering
sc = yt.create_scene(ds, field=('gas', 'density'))

# Now increase the resolution
sc.camera.resolution = (1024, 1024)

# Set the camera focus to a position that is offset from the center of
# the domain
sc.camera.focus = ds.arr([0.3, 0.3, 0.3], 'unitary')

# Move the camera position to the other side of the dataset
sc.camera.position = ds.arr([0, 0, 0], 'unitary')

# save to disk with a custom filename and apply sigma clipping to eliminate
# very bright pixels, producing an image with better contrast.
sc.render()
sc.save('custom.png', sigma_clip=4)
../_images/custom_camera_volume_rendering__custom.png

Volume Rendering with a Custom Transfer Function

In this recipe we modify the Simple Volume Rendering recipe to use customized camera properties. See 3D Visualization and Volume Rendering for more information.

(custom_transfer_function_volume_rendering.py)

import yt
import numpy as np

# Load the dataset
ds = yt.load("Enzo_64/DD0043/data0043")

# Create a volume rendering
sc = yt.create_scene(ds, field=('gas', 'density'))

# Modify the transfer function

# First get the render source, in this case the entire domain, with field ('gas','density')
render_source = sc.get_source(0)

# Clear the transfer function
render_source.transfer_function.clear()

# Map a range of density values (in log space) to the Reds_r colormap
render_source.transfer_function.map_to_colormap(
    np.log10(ds.quan(5.0e-31, 'g/cm**3')),
    np.log10(ds.quan(1.0e-29, 'g/cm**3')),
    scale=30.0, colormap='RdBu_r')

sc.save('new_tf.png')
../_images/custom_transfer_function_volume_rendering__new_tf.png

Volume Rendering with Sigma Clipping

In this recipe we output several images with different values of sigma_clip set in order to change the contrast of the resulting image. See Improving Image Contrast with Sigma Clipping for more information.

(sigma_clip.py)

import yt

# Load the dataset.
ds = yt.load("enzo_tiny_cosmology/RD0009/RD0009")

# Create a volume rendering, which will determine data bounds, use the first
# acceptable field in the field_list, and set up a default transfer function.

# Render and save output images with different levels of sigma clipping.
# Sigma clipping removes the highest intensity pixels in a volume render,
# which affects the overall contrast of the image.
sc = yt.create_scene(ds, field=('gas', 'density'))
sc.render()
sc.save('clip_0.png')
sc.save('clip_2.png', sigma_clip=2)
sc.save('clip_4.png', sigma_clip=4)
sc.save('clip_6.png', sigma_clip=6)
../_images/sigma_clip__clip_0.png ../_images/sigma_clip__clip_2.png ../_images/sigma_clip__clip_4.png ../_images/sigma_clip__clip_6.png

Zooming into an Image

This is a recipe that takes a slice through the most dense point, then creates a bunch of frames as it zooms in. It’s important to note that this particular recipe is provided to show how to be more flexible and add annotations and the like – the base system, of a zoomin, is provided by the “yt zoomin” command on the command line. See Slice Plots and Plot Modifications: Overplotting Contours, Velocities, Particles, and More for more information.

(zoomin_frames.py)

import yt
import numpy as np

# load data
fn = "IsolatedGalaxy/galaxy0030/galaxy0030"
ds = yt.load(fn)

# This is the number of frames to make -- below, you can see how this is used.
n_frames = 5

# This is the minimum size in smallest_dx of our last frame.
# Usually it should be set to something like 400, but for THIS
# dataset, we actually don't have that great of resolution.
min_dx = 40

frame_template = "frame_%05i" # Template for frame filenames

p = yt.SlicePlot(ds, "z", "density") # Add our slice, along z
p.annotate_contour("temperature") # We'll contour in temperature

# What we do now is a bit fun.  "enumerate" returns a tuple for every item --
# the index of the item, and the item itself.  This saves us having to write
# something like "i = 0" and then inside the loop "i += 1" for ever loop.  The
# argument to enumerate is the 'logspace' function, which takes a minimum and a
# maximum and the number of items to generate.  It returns 10^power of each
# item it generates.

for i,v in enumerate(np.logspace(0,
                                 np.log10(ds.index.get_smallest_dx()*min_dx),
                                 n_frames)):
    # We set our width as necessary for this frame
    p.set_width(v, 'unitary')
    # save
    p.save(frame_template % (i))
../_images/zoomin_frames__frame_00000_Slice_z_density.png ../_images/zoomin_frames__frame_00001_Slice_z_density.png ../_images/zoomin_frames__frame_00002_Slice_z_density.png ../_images/zoomin_frames__frame_00003_Slice_z_density.png ../_images/zoomin_frames__frame_00004_Slice_z_density.png

Various Lens Types for Volume Rendering

This example illustrates the usage and feature of different lenses for volume rendering.

(various_lens.py)

import yt
from yt.visualization.volume_rendering.api import Scene, VolumeSource
import numpy as np

field = ("gas", "density")

# normal_vector points from camera to the center of tbe final projection.
# Now we look at the positive x direction.
normal_vector = [1., 0., 0.]
# north_vector defines the "top" direction of the projection, which is
# positive z direction here.
north_vector = [0., 0., 1.]

# Follow the simple_volume_rendering cookbook for the first part of this.
ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
sc = Scene()
vol = VolumeSource(ds, field=field)
tf = vol.transfer_function
tf.grey_opacity = True

# Plane-parallel lens
cam = sc.add_camera(ds, lens_type='plane-parallel')
# Set the resolution of tbe final projection.
cam.resolution = [250, 250]
# Set the location of the camera to be (x=0.2, y=0.5, z=0.5)
# For plane-parallel lens, the location info along the normal_vector (here
# is x=0.2) is ignored.
cam.position = ds.arr(np.array([0.2, 0.5, 0.5]), 'code_length')
# Set the orientation of the camera.
cam.switch_orientation(normal_vector=normal_vector,
                       north_vector=north_vector)
# Set the width of the camera, where width[0] and width[1] specify the length and
# height of final projection, while width[2] in plane-parallel lens is not used.
cam.set_width(ds.domain_width * 0.5)
sc.add_source(vol)
sc.render()
sc.save('lens_plane-parallel.png', sigma_clip=6.0)

# Perspective lens
cam = sc.add_camera(ds, lens_type='perspective')
cam.resolution = [250, 250]
# Standing at (x=0.2, y=0.5, z=0.5), we look at the area of x>0.2 (with some open angle
# specified by camera width) along the positive x direction.
cam.position = ds.arr([0.2, 0.5, 0.5], 'code_length')
cam.switch_orientation(normal_vector=normal_vector,
                       north_vector=north_vector)
# Set the width of the camera, where width[0] and width[1] specify the length and
# height of the final projection, while width[2] specifies the distance between the
# camera and the final image.
cam.set_width(ds.domain_width * 0.5)
sc.add_source(vol)
sc.render()
sc.save('lens_perspective.png', sigma_clip=6.0)

# Stereo-perspective lens
cam = sc.add_camera(ds, lens_type='stereo-perspective')
# Set the size ratio of the final projection to be 2:1, since stereo-perspective lens
# will generate the final image with both left-eye and right-eye ones jointed together.
cam.resolution = [500, 250]
cam.position = ds.arr([0.2, 0.5, 0.5], 'code_length')
cam.switch_orientation(normal_vector=normal_vector,
                       north_vector=north_vector)
cam.set_width(ds.domain_width*0.5)
# Set the distance between left-eye and right-eye.
cam.lens.disparity = ds.domain_width[0] * 1.e-3
sc.add_source(vol)
sc.render()
sc.save('lens_stereo-perspective.png', sigma_clip=6.0)

# Fisheye lens
dd = ds.sphere(ds.domain_center, ds.domain_width[0] / 10)
cam = sc.add_camera(dd, lens_type='fisheye')
cam.resolution = [250, 250]
v, c = ds.find_max(field)
cam.set_position(c - 0.0005 * ds.domain_width)
cam.switch_orientation(normal_vector=normal_vector,
                       north_vector=north_vector)
cam.set_width(ds.domain_width)
cam.lens.fov = 360.0
sc.add_source(vol)
sc.render()
sc.save('lens_fisheye.png', sigma_clip=6.0)

# Spherical lens
cam = sc.add_camera(ds, lens_type='spherical')
# Set the size ratio of the final projection to be 2:1, since spherical lens
# will generate the final image with length of 2*pi and height of pi.
# Recommended resolution for YouTube 360-degree videos is [3840, 2160]
cam.resolution = [500, 250]
# Standing at (x=0.4, y=0.5, z=0.5), we look in all the radial directions
# from this point in spherical coordinate.
cam.position = ds.arr([0.4, 0.5, 0.5], 'code_length')
cam.switch_orientation(normal_vector=normal_vector,
                       north_vector=north_vector)
# In (stereo)spherical camera, camera width is not used since the entire volume
# will be rendered
sc.add_source(vol)
sc.render()
sc.save('lens_spherical.png', sigma_clip=6.0)

# Stereo-spherical lens
cam = sc.add_camera(ds, lens_type='stereo-spherical')
# Set the size ratio of the final projection to be 1:1, since spherical-perspective lens
# will generate the final image with both left-eye and right-eye ones jointed together,
# with left-eye image on top and right-eye image on bottom.
# Recommended resolution for YouTube virtual reality videos is [3840, 2160]
cam.resolution = [500, 500]
cam.position = ds.arr([0.4, 0.5, 0.5], 'code_length')
cam.switch_orientation(normal_vector=normal_vector,
                       north_vector=north_vector)
# In (stereo)spherical camera, camera width is not used since the entire volume
# will be rendered
# Set the distance between left-eye and right-eye.
cam.lens.disparity = ds.domain_width[0] * 1.e-3
sc.add_source(vol)
sc.render()
sc.save('lens_stereo-spherical.png', sigma_clip=6.0)
../_images/various_lens__lens_fisheye.png ../_images/various_lens__lens_perspective.png ../_images/various_lens__lens_plane-parallel.png ../_images/various_lens__lens_spherical.png ../_images/various_lens__lens_stereo-perspective.png ../_images/various_lens__lens_stereo-spherical.png

Opaque Volume Rendering

This recipe demonstrates how to make semi-opaque volume renderings, but also how to step through and try different things to identify the type of volume rendering you want. See Opacity for more information.

(opaque_rendering.py)

import yt
import numpy as np

ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")

# We start by building a default volume rendering scene

im, sc = yt.volume_render(ds, field=("gas","density"), fname="v0.png", sigma_clip=6.0)

sc.camera.set_width(ds.arr(0.1,'code_length'))
tf = sc.get_source(0).transfer_function
tf.clear()
tf.add_layers(4, 0.01, col_bounds = [-27.5,-25.5],
        alpha=np.logspace(-3,0,4), colormap = 'RdBu_r')
sc.render()
sc.save("v1.png", sigma_clip=6.0)

# In this case, the default alphas used (np.logspace(-3,0,Nbins)) does not
# accentuate the outer regions of the galaxy. Let's start by bringing up the
# alpha values for each contour to go between 0.1 and 1.0

tf = sc.get_source(0).transfer_function
tf.clear()
tf.add_layers(4, 0.01, col_bounds = [-27.5,-25.5],
        alpha=np.logspace(0,0,4), colormap = 'RdBu_r')
sc.render()
sc.save("v2.png", sigma_clip=6.0)

# Now let's set the grey_opacity to True.  This should make the inner portions
# start to be obcured

tf.grey_opacity = True
sc.render()
sc.save("v3.png", sigma_clip=6.0)

# That looks pretty good, but let's start bumping up the opacity.

tf.clear()
tf.add_layers(4, 0.01, col_bounds = [-27.5,-25.5],
        alpha=10.0*np.ones(4,dtype='float64'), colormap = 'RdBu_r')
sc.render()
sc.save("v4.png", sigma_clip=6.0)

# Let's bump up again to see if we can obscure the inner contour.

tf.clear()
tf.add_layers(4, 0.01, col_bounds = [-27.5,-25.5],
        alpha=30.0*np.ones(4,dtype='float64'), colormap = 'RdBu_r')
sc.render()
sc.save("v5.png", sigma_clip=6.0)

# Now we are losing sight of everything.  Let's see if we can obscure the next
# layer

tf.clear()
tf.add_layers(4, 0.01, col_bounds = [-27.5,-25.5],
        alpha=100.0*np.ones(4,dtype='float64'), colormap = 'RdBu_r')
sc.render()
sc.save("v6.png", sigma_clip=6.0)

# That is very opaque!  Now lets go back and see what it would look like with
# grey_opacity = False

tf.grey_opacity=False
sc.render()
sc.save("v7.png", sigma_clip=6.0)

# That looks pretty different, but the main thing is that you can see that the
# inner contours are somewhat visible again.
../_images/opaque_rendering__v0.png ../_images/opaque_rendering__v1.png ../_images/opaque_rendering__v2.png ../_images/opaque_rendering__v3.png ../_images/opaque_rendering__v4.png ../_images/opaque_rendering__v5.png ../_images/opaque_rendering__v6.png ../_images/opaque_rendering__v7.png

Downsampling Data for Volume Rendering

This recipe demonstrates how to downsample data in a simulation to speed up volume rendering. See 3D Visualization and Volume Rendering for more information.

(amrkdtree_downsampling.py)

# Using AMRKDTree Homogenized Volumes to examine large datasets
# at lower resolution.

# In this example we will show how to use the AMRKDTree to take a simulation
# with 8 levels of refinement and only use levels 0-3 to render the dataset.

# Currently this cookbook is flawed in that the data that is covered by the
# higher resolution data gets masked during the rendering.  This should be
# fixed by changing either the data source or the code in
# yt/utilities/amr_kdtree.py where data is being masked for the partitioned
# grid.  Right now the quick fix is to create a data_collection, but this
# will only work for patch based simulations that have ds.index.grids.

# We begin by loading up yt, and importing the AMRKDTree
import numpy as np

import yt
from yt.utilities.amr_kdtree.api import AMRKDTree

# Load up a dataset and define the kdtree
ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
im, sc = yt.volume_render(ds, 'density', fname='v0.png')
sc.camera.set_width(ds.arr(100, 'kpc'))
render_source = sc.get_source(0)
kd=render_source.volume

# Print out specifics of KD Tree
print("Total volume of all bricks = %i" % kd.count_volume())
print("Total number of cells = %i" % kd.count_cells())

new_source = ds.all_data()
new_source.max_level=3
kd_low_res = AMRKDTree(ds, data_source=new_source)
print(kd_low_res.count_volume())
print(kd_low_res.count_cells())

# Now we pass this in as the volume to our camera, and render the snapshot
# again.

render_source.set_volume(kd_low_res)
render_source.set_field('density')
sc.render()
sc.save("v1.png", sigma_clip=6.0)

# This operation was substantiall faster.  Now lets modify the low resolution
# rendering until we find something we like.

tf = render_source.transfer_function
tf.clear()
tf.add_layers(4, 0.01, col_bounds=[-27.5, -25.5],
              alpha=np.ones(4, dtype='float64'), colormap='RdBu_r')
sc.render()
sc.save("v2.png", sigma_clip=6.0)

# This looks better.  Now let's try turning on opacity.

tf.grey_opacity = True
sc.render()
sc.save("v3.png", sigma_clip=6.0)
#
## That seemed to pick out som interesting structures.  Now let's bump up the
## opacity.
#
tf.clear()
tf.add_layers(4, 0.01, col_bounds=[-27.5, -25.5],
              alpha=10.0 * np.ones(4, dtype='float64'), colormap='RdBu_r')
sc.render()
sc.save("v4.png", sigma_clip=6.0)
#
## This looks pretty good, now lets go back to the full resolution AMRKDTree
#
render_source.set_volume(kd)
sc.render()
sc.save("v5.png", sigma_clip=6.0)

# This looks great!
../_images/amrkdtree_downsampling__v0.png ../_images/amrkdtree_downsampling__v1.png ../_images/amrkdtree_downsampling__v2.png ../_images/amrkdtree_downsampling__v3.png ../_images/amrkdtree_downsampling__v4.png ../_images/amrkdtree_downsampling__v5.png

Volume Rendering with Bounding Box and Overlaid Grids

This recipe demonstrates how to overplot a bounding box on a volume rendering as well as overplotting grids representing the level of refinement achieved in different regions of the code. See Annotations for more information.

(rendering_with_box_and_grids.py)

import yt

# Load the dataset.
ds = yt.load("Enzo_64/DD0043/data0043")
sc = yt.create_scene(ds, ('gas', 'density'))

# You may need to adjust the alpha values to get a rendering with good contrast
# For annotate_domain, the fourth color value is alpha.

# Draw the domain boundary
sc.annotate_domain(ds, color=[1, 1, 1, 0.01])
sc.save("%s_vr_domain.png" % ds, sigma_clip=4)

# Draw the grid boundaries
sc.annotate_grids(ds, alpha=0.01)
sc.save("%s_vr_grids.png" % ds, sigma_clip=4)

# Draw a coordinate axes triad
sc.annotate_axes(alpha=0.01)
sc.save("%s_vr_coords.png" % ds, sigma_clip=4)
../_images/rendering_with_box_and_grids__data0043_vr_coords.png ../_images/rendering_with_box_and_grids__data0043_vr_domain.png ../_images/rendering_with_box_and_grids__data0043_vr_grids.png

Volume Rendering with Annotation

This recipe demonstrates how to write the simulation time, show an axis triad indicating the direction of the coordinate system, and show the transfer function on a volume rendering. Please note that this recipe relies on the old volume rendering interface. While one can continue to use this interface, it may be incompatible with some of the new developments and the infrastructure described in 3D Visualization and Volume Rendering.

(vol-annotated.py)

import yt

ds = yt.load('Enzo_64/DD0043/data0043')

sc = yt.create_scene(ds, lens_type='perspective')

source = sc[0]

source.set_field('density')
source.set_log(True)

# Set up the camera parameters: focus, width, resolution, and image orientation
sc.camera.focus = ds.domain_center
sc.camera.resolution = 1024
sc.camera.north_vector = [0, 0, 1]
sc.camera.position = [1.7, 1.7, 1.7]

# You may need to adjust the alpha values to get an image with good contrast.
# For the annotate_domain call, the fourth value in the color tuple is the
# alpha value.
sc.annotate_axes(alpha=.02)
sc.annotate_domain(ds, color=[1, 1, 1, .01])

text_string = "T = {} Gyr".format(float(ds.current_time.to('Gyr')))

# save an annotated version of the volume rendering including a representation
# of the transfer function and a nice label showing the simulation time.
sc.save_annotated("vol_annotated.png", sigma_clip=6,
                  text_annotate=[[(.1, 1.05), text_string]])
../_images/vol-annotated__vol_annotated.png

Volume Rendering with Points

This recipe demonstrates how to make a volume rendering composited with point sources. This could represent star or dark matter particles, for example.

(vol-points.py)

import yt
import numpy as np
from yt.visualization.volume_rendering.api import PointSource
from yt.units import kpc

ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')

sc = yt.create_scene(ds)

np.random.seed(1234567)

npoints = 1000

# Random particle positions
vertices = np.random.random([npoints, 3])*200*kpc

# Random colors
colors = np.random.random([npoints, 4])

# Set alpha value to something that produces a good contrast with the volume
# rendering
colors[:, 3] = 0.1

points = PointSource(vertices, colors=colors)
sc.add_source(points)

sc.camera.width = 300*kpc

sc.save(sigma_clip=5)
../_images/vol-points__galaxy0030_Render_density.png

Volume Rendering with Lines

This recipe demonstrates how to make a volume rendering composited with line sources.

(vol-lines.py)

import yt
import numpy as np
from yt.visualization.volume_rendering.api import LineSource
from yt.units import kpc

ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')

sc = yt.create_scene(ds)

np.random.seed(1234567)

nlines = 50
vertices = (np.random.random([nlines, 2, 3]) - 0.5) * 200 * kpc
colors = np.random.random([nlines, 4])
colors[:, 3] = 0.1

lines = LineSource(vertices, colors)
sc.add_source(lines)

sc.camera.width = 300*kpc

sc.save(sigma_clip=4.0)
../_images/vol-lines__galaxy0030_Render_density.png

Advanced Interactive Data Visualization

This recipe demonstrates how to manually create all components required to start the Interactive Data Visualization. For more information see Interactive Data Visualization.

(opengl_vr.py)

import yt
from yt.visualization.volume_rendering.interactive_vr import \
    SceneGraph, BlockCollection, TrackballCamera
from yt.visualization.volume_rendering.interactive_loop import \
    RenderingContext

ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")

# Create GLUT window
rc = RenderingContext(1280, 960)

# Create a 3d Texture from all_data()
collection = BlockCollection()
dd = ds.all_data()
collection.add_data(dd, "density")

# Initiliaze basic Scene and pass the data
scene = SceneGraph()
scene.add_collection(collection)

# Create default camera
position = (1.0, 1.0, 1.0)
c = TrackballCamera(position=position, focus=ds.domain_center,
                    near_plane=0.1)

# Start rendering loop
rc.start_loop(scene, c)

Embedding Interactive Data Visualization

This recipe demonstrates how to embed the Interactive Data Visualization inside the Jupyter notebook. For more information see Interactive Data Visualization.

(opengl_ipython.py)

import yt
from yt.visualization.volume_rendering.interactive_vr import \
    SceneGraph, BlockCollection, TrackballCamera
from yt.visualization.volume_rendering.interactive_loop import \
    RenderingContext
from yt.visualization.volume_rendering import glfw_inputhook 

rc = RenderingContext(1280, 960)

scene = SceneGraph()
collection = BlockCollection()

ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")

dd = ds.all_data()
collection.add_data(dd, "density")

scene.add_collection(collection)

position = (1.0, 1.0, 1.0)
c = TrackballCamera(position=position, focus=ds.domain_center,
                    near_plane=0.1)

callbacks = rc.setup_loop(scene, c)
rl = rc(scene, c, callbacks)

# To make this work from IPython execute:
#
# glfw_inputhook.inputhook_manager.enable_gui("glfw", app=rl)

Plotting Streamlines

This recipe demonstrates how to display streamlines in a simulation. (Note: streamlines can also be queried for values!) See Streamlines: Tracking the Trajectories of Tracers in your Data for more information.

(streamlines.py)

import yt
import numpy as np
import matplotlib.pylab as pl

from yt.visualization.api import Streamlines
from yt.units import Mpc
from mpl_toolkits.mplot3d import Axes3D

# Load the dataset
ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')

# Define c: the center of the box, N: the number of streamlines,
# scale: the spatial scale of the streamlines relative to the boxsize,
# and then pos: the random positions of the streamlines.
c = ds.domain_center
N = 100
scale = ds.domain_width[0]
pos_dx = np.random.random((N,3))*scale-scale/2.
pos = c+pos_dx

# Create streamlines of the 3D vector velocity and integrate them through
# the box defined above
streamlines = Streamlines(ds, pos, 'velocity_x', 'velocity_y', 'velocity_z',
                          length=1.0*Mpc, get_magnitude=True)
streamlines.integrate_through_volume()

# Create a 3D plot, trace the streamlines throught the 3D volume of the plot
fig=pl.figure()
ax = Axes3D(fig)
for stream in streamlines.streamlines:
    stream = stream[np.all(stream != 0.0, axis=1)]
    ax.plot3D(stream[:,0], stream[:,1], stream[:,2], alpha=0.1)

# Save the plot to disk.
pl.savefig('streamlines.png')
../_images/streamlines__streamlines.png

Plotting Isocontours

This recipe demonstrates how to extract an isocontour and then plot it in matplotlib, coloring the surface by a second quantity. See 3D Surfaces and Sketchfab for more information.

(surface_plot.py)

import yt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import matplotlib.pyplot as plt

# Load the dataset
ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")

# Create a sphere object centered on the highest density point in the simulation
# with radius 1 Mpc
sphere = ds.sphere("max", (1.0, "Mpc"))

# Identify the isodensity surface in this sphere with density = 1e-24 g/cm^3
surface = ds.surface(sphere, "density", 1e-24)

# Color this isodensity surface according to the log of the temperature field
colors = yt.apply_colormap(np.log10(surface["temperature"]), cmap_name="hot")

# Create a 3D matplotlib figure for visualizing the surface
fig = plt.figure()
ax = fig.gca(projection='3d')
p3dc = Poly3DCollection(surface.triangles, linewidth=0.0)

# Set the surface colors in the right scaling [0,1]
p3dc.set_facecolors(colors[0,:,:]/255.)
ax.add_collection(p3dc)

# Let's keep the axis ratio fixed in all directions by taking the maximum
# extent in one dimension and make it the bounds in all dimensions
max_extent = (surface.vertices.max(axis=1) - surface.vertices.min(axis=1)).max()
centers = (surface.vertices.max(axis=1) + surface.vertices.min(axis=1)) / 2
bounds = np.zeros([3,2])
bounds[:,0] = centers[:] - max_extent/2
bounds[:,1] = centers[:] + max_extent/2
ax.auto_scale_xyz(bounds[0,:], bounds[1,:], bounds[2,:])

# Save the figure
plt.savefig("%s_Surface.png" % ds)
../_images/surface_plot__galaxy0030_Surface.png

Plotting Isocontours and Streamlines

This recipe plots both isocontours and streamlines simultaneously. Note that this will not include any blending, so streamlines that are occluded by the surface will still be visible. See Streamlines: Tracking the Trajectories of Tracers in your Data and 3D Surfaces and Sketchfab for more information.

(streamlines_isocontour.py)

import yt
import numpy as np
from yt.visualization.api import Streamlines
import matplotlib.pylab as pl
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

# Load the dataset
ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")

# Define c: the center of the box, N: the number of streamlines,
# scale: the spatial scale of the streamlines relative to the boxsize,
# and then pos: the random positions of the streamlines.
c = ds.arr([0.5]*3, 'code_length')
N = 30
scale = ds.quan(15, 'kpc').in_units('code_length')  # 15 kpc in code units
pos_dx = np.random.random((N,3))*scale-scale/2.
pos = c+pos_dx

# Create the streamlines from these positions with the velocity fields as the
# fields to be traced
streamlines = Streamlines(ds, pos, 'velocity_x', 'velocity_y', 'velocity_z', length=1.0)
streamlines.integrate_through_volume()

# Create a 3D matplotlib figure for visualizing the streamlines
fig=pl.figure()
ax = Axes3D(fig)

# Trace the streamlines through the volume of the 3D figure
for stream in streamlines.streamlines:
    stream = stream[np.all(stream != 0.0, axis=1)]

    # Make the colors of each stream vary continuously from blue to red
    # from low-x to high-x of the stream start position (each color is R, G, B)
    # can omit and just set streamline colors to a fixed color
    x_start_pos = ds.arr(stream[0,0], 'code_length')
    x_start_pos -= ds.arr(0.5, 'code_length')
    x_start_pos /= scale
    x_start_pos += 0.5
    color = np.array([x_start_pos, 0, 1-x_start_pos])

    # Plot the stream in 3D
    ax.plot3D(stream[:,0], stream[:,1], stream[:,2], alpha=0.3, color=color)

# Create a sphere object centered on the highest density point in the simulation
# with radius = 1 Mpc
sphere = ds.sphere("max", (1.0, "Mpc"))

# Identify the isodensity surface in this sphere with density = 1e-24 g/cm^3
surface = ds.surface(sphere, "density", 1e-24)

# Color this isodensity surface according to the log of the temperature field
colors = yt.apply_colormap(np.log10(surface["temperature"]), cmap_name="hot")

# Render this surface
p3dc = Poly3DCollection(surface.triangles, linewidth=0.0)
colors = colors[0,:,:]/255.  # scale to [0,1]
colors[:,3] = 0.3            # alpha = 0.3
p3dc.set_facecolors(colors)
ax.add_collection(p3dc)

# Save the figure
pl.savefig('streamlines_isocontour.png')
../_images/streamlines_isocontour__streamlines_isocontour.png