## Saturday, August 1, 2015

### Illuminating Surface Plots

Matplotlib provides functions for visualizing three-dimensional data sets. One useful tool is a surface plot. A surface plot is a two-dimensional projection of a three-dimensional object. Much like a sketch artist, Python uses techniques like perspective and shading to give the illusion of a three-dimensional object in space. In this post, I describe how you can control the lighting of a surface plot.

## Surface Plots

First, let’s look at some of the options available with the default three-dimensional plotting tools. This script will create a surface plot of a Bessel function. Its ripples will emphasize the effects of lighting later.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D     # Import 3D plotting tools.
from scipy.special import jn                # Import Bessel function.

# Define grid of points.
points = np.linspace(-10, 10, 51)
X, Y = np.meshgrid(points, points)
R = np.sqrt(X**2 + Y**2)
Z = jn(0,R)

# Create 3D surface plot.
ax = Axes3D(plt.figure())
ax.plot_surface(X, Y, Z, rstride=1, cstride=1)


The default surface plot is a single color with grid lines and some shading. To get rid of the grid lines use the following plot commands instead:

ax = Axes3D(plt.figure())
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0, antialiased=False,
color='green')


We can use a colormap to assign colors to the figure based on the height of the surface

ax = Axes3D(plt.figure())
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0, antialiased=False,
cmap='coolwarm')


Assigning colors according to height may not be what you want, and there is no shading when using a color map. Furthermore, when you specify a single color, you cannot adjust the lighting angle to produce different shading effects. Sometimes, you may want to control these lighting effects, and Matplotlib provides a way.

## Turn on the Lights

Tucked away in Matplotlib is an object called LightSource. It allows you to simulate illuminating a surface using a virtual light source placed at a location of your choosing. (LightSource creates an “illuminated intensity map.” You can find the details of the model in the source code.) It does not provide the same control or as many features as the lighting tools of commercial packages like MATLAB or Mathematica, but it is sufficient to produce some nice plots.

To use the object, import it from matplotlib.colors and then apply its shade method to the data set:

# Get lighting object for shading surface plots.
from matplotlib.colors import LightSource

# Get colormaps to use with lighting object.
from matplotlib import cm

# Create an instance of a LightSource and use it to illuminate the surface.
light = LightSource(90, 45)


The two arguments to LightSource are the azimuth and elevation angles of the light source. (0,0) corresponds to a light placed along the x-axis. As the name implies, the elevation is the angle above the xy-plane in degrees. The virtual light source position is then rotated about the vertical axis by the azimuth angle (also in degrees). (Don’t confuse these parameters with the similarly named parameters specifying the observer’s position!)

The function requires a single argument: a two-dimensional array — here, Z. LightSource interprets each data point as the height of a surface above a point in the xy-plane. It also assumes these points have the same spacing in the x and y directions. If you are not using Cartesian coordinates and uniform spacing, you may be surprised by the result.

The object returned by light.shade is a NumPy array of RGBA values for each data point. (For each point in the input array, light.shade returns a 4-element array of the Red, Green, Blue, and Alpha value for that point. Alpha controls the transparency of the point.) Other plotting tools can use this data to draw a shaded surface. To use this array instead of the color or cmap options of the surface plotting command, pass the array with a keyword argument:

ax = Axes3D(plt.figure())
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0, antialiased=False,
facecolors=illuminated_surface)


If you prefer to shade a surface of uniform color instead of using a color map, or if you have a colored surface that you wish to shade, a LightSource object offers a second method called shade_rgb. You have to pass the function two arguments: an array of RGB values and a data set giving the height of each point.

As an example, let’s transform a white surface so we can see the shading effects independent of any coloring. The RGB values for white are [1.0, 1.0, 1.0]. (Red, green, and blue values are all at maximum.) To create a uniform white surface, we need to create an array with three elements for every element of the data set Z, with each entry set to 1.0. The following code will create the RGB array, shade it, and plot it:

rgb = np.ones((Z.shape[0], Z.shape[1], 3))
ax = Axes3D(plt.figure())
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0, antialiased=False,
facecolors=illuminated_surface)


To change the color of the shaded surface, we can use NumPy array math. Just make a three-element array with the RGB values of your target color and multiply rgb by this array before shading.

# Create a shaded green surface.
green = np.array([0,1.0,0])
green_surface = light.shade_rgb(rgb * green, Z)
ax = Axes3D(plt.figure())
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0, antialiased=False,
facecolors=green_surface)


The figure below illustrates the techniques described here. The same surface is shown with four different color and lighting configurations. The code that produced the figures is also included below.

# =========================================================================
# Author:   Jesse M. Kinder
# Created:  2015 Jul 27
# Modified: 2015 Jul 31
# -------------------------------------------------------------------------
# Demonstrate shading of surface plots using Matplotlib's LightSource.
# -------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Import Bessel function.
from scipy.special import jn

# Import colormaps.
from matplotlib import cm

# Import lighting object for shading surface plots.
from matplotlib.colors import LightSource

# Define grid of points.
points = np.linspace(-10, 10, 101)
X, Y = np.meshgrid(points, points)
R = np.sqrt(X**2 + Y**2)
Z = jn(0,R)

# Create an rgb array for single-color surfaces.
white = np.ones((Z.shape[0], Z.shape[1], 3))
red = white * np.array([1,0,0])
green = white * np.array([0,1,0])
blue = white * np.array([0,0,1])

# Set view parameters for all subplots.
azimuth = 45
altitude = 60

# Create empty figure.
fig = plt.figure(figsize=(18,12))

# -------------------------------------------------------------------------
# Generate first subplot.
# -------------------------------------------------------------------------
# Create a light source object for light from
# 0 degrees azimuth, 0 degrees elevation.
light = LightSource(0, 0)

# Generate face colors for a shaded surface using either
# a color map or the uniform rgb color specified above.

# Create a subplot with 3d plotting capabilities.
# This command will fail if Axes3D was not imported.
ax.view_init(altitude, azimuth)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0,
antialiased=False, facecolors=illuminated_surface)

# -------------------------------------------------------------------------
# Repeat the commands above for the other three subplots, but use different
# illumination angles and colors.
# -------------------------------------------------------------------------
light = LightSource(90, 0)

ax.view_init(altitude, azimuth)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0,
antialiased=False, facecolors=illuminated_surface)

# -------------------------------------------------------------------------
light = LightSource(90, 45)

ax.view_init(altitude, azimuth)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0,
antialiased=False, facecolors=illuminated_surface)

# -------------------------------------------------------------------------
light = LightSource(180, 45)

ax.view_init(altitude, azimuth)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0,
antialiased=False, facecolors=illuminated_surface)

# -------------------------------------------------------------------------
plt.tight_layout()