matplotlib · jenshnielsen · Sep 29, 2014 · Jul 21, 2014 · Jul 21, 2014 · Jul 21, 2014
diff --git a/CHANGELOG b/CHANGELOG
@@ -1,3 +1,9 @@
+2014-09-27 Overhauled `colors.LightSource`.  Added `LightSource.hillshade` to
+	   allow the independent generation of illumination maps. Added new
+	   types of blending for creating more visually appealing shaded relief
+	   plots (e.g.  `blend_mode="overlay"`, etc, in addition to the legacy
+	   "hsv" mode).
+
 2014-06-10 Added Colorbar.remove()
 
 2014-06-07 Fixed bug so radial plots can be saved as ps in py3k.

diff --git a/examples/mplot3d/custom_shaded_3d_surface.py b/examples/mplot3d/custom_shaded_3d_surface.py
@@ -0,0 +1,32 @@
+"""
+Demonstrates using custom hillshading in a 3D surface plot.
+"""
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib import cbook
+from matplotlib import cm
+from matplotlib.colors import LightSource
+import matplotlib.pyplot as plt
+import numpy as np
+
+filename = cbook.get_sample_data('jacksboro_fault_dem.npz', asfileobj=False)
+with np.load(filename) as dem:
+    z = dem['elevation']
+    nrows, ncols = z.shape
+    x = np.linspace(dem['xmin'], dem['xmax'], ncols)
+    y = np.linspace(dem['ymin'], dem['ymax'], nrows)
+    x, y = np.meshgrid(x, y)
+
+region = np.s_[5:50, 5:50]
+x, y, z = x[region], y[region], z[region]
+
+fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))
+
+ls = LightSource(270, 45)
+# To use a custom hillshading mode, override the built-in shading and pass
+# in the rgb colors of the shaded surface calculated from "shade".
+rgb = ls.shade(z, cmap=cm.gist_earth, vert_exag=0.1, blend_mode='soft')
+surf = ax.plot_surface(x, y, z, rstride=1, cstride=1, facecolors=rgb,
+                       linewidth=0, antialiased=False, shade=False)
+
+plt.show()
+
diff --git a/examples/pylab_examples/shading_example.py b/examples/pylab_examples/shading_example.py
@@ -1,28 +1,55 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib.colors import LightSource
+from matplotlib.cbook import get_sample_data
 
-# example showing how to make shaded relief plots 
+# Example showing how to make shaded relief plots
 # like Mathematica
 # (http://reference.wolfram.com/mathematica/ref/ReliefPlot.html)
 # or Generic Mapping Tools
 # (http://gmt.soest.hawaii.edu/gmt/doc/gmt/html/GMT_Docs/node145.html)
 
-# test data
-X,Y=np.mgrid[-5:5:0.05,-5:5:0.05]
-Z=np.sqrt(X**2+Y**2)+np.sin(X**2+Y**2)
-# create light source object.
-ls = LightSource(azdeg=0,altdeg=65)
-# shade data, creating an rgb array.
-rgb = ls.shade(Z,plt.cm.copper)
-# plot un-shaded and shaded images.
-plt.figure(figsize=(12,5))
-plt.subplot(121)
-plt.imshow(Z,cmap=plt.cm.copper)
-plt.title('imshow')
-plt.xticks([]); plt.yticks([])
-plt.subplot(122)
-plt.imshow(rgb)
-plt.title('imshow with shading')
-plt.xticks([]); plt.yticks([])
-plt.show()
+def main():
+    # Test data
+    x, y = np.mgrid[-5:5:0.05, -5:5:0.05]
+    z = 5 * (np.sqrt(x**2 + y**2) + np.sin(x**2 + y**2))
+
+    filename = get_sample_data('jacksboro_fault_dem.npz', asfileobj=False)
+    with np.load(filename) as dem:
+        elev = dem['elevation']
+
+    fig = compare(z, plt.cm.copper)
+    fig.suptitle('HSV Blending Looks Best with Smooth Surfaces', y=0.95)
+
+    fig = compare(elev, plt.cm.gist_earth, ve=0.05)
+    fig.suptitle('Overlay Blending Looks Best with Rough Surfaces', y=0.95)
+
+    plt.show()
+
+def compare(z, cmap, ve=1):
+    # Create subplots and hide ticks
+    fig, axes = plt.subplots(ncols=2, nrows=2)
+    for ax in axes.flat:
+        ax.set(xticks=[], yticks=[])
+
+    # Illuminate the scene from the northwest
+    ls = LightSource(azdeg=315, altdeg=45)
+
+    axes[0, 0].imshow(z, cmap=cmap)
+    axes[0, 0].set(xlabel='Colormapped Data')
+
+    axes[0, 1].imshow(ls.hillshade(z, vert_exag=ve), cmap='gray')
+    axes[0, 1].set(xlabel='Illumination Intensity')
+
+    rgb = ls.shade(z, cmap=cmap, vert_exag=ve, blend_mode='hsv')
+    axes[1, 0].imshow(rgb)
+    axes[1, 0].set(xlabel='Blend Mode: "hsv" (default)')
+
+    rgb = ls.shade(z, cmap=cmap, vert_exag=ve, blend_mode='overlay')
+    axes[1, 1].imshow(rgb)
+    axes[1, 1].set(xlabel='Blend Mode: "overlay"')
+
+    return fig
+
+if __name__ == '__main__':
+    main()
diff --git a/examples/specialty_plots/advanced_hillshading.py b/examples/specialty_plots/advanced_hillshading.py
@@ -0,0 +1,68 @@
+"""
+Demonstrates a few common tricks with shaded plots.
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.colors import LightSource, Normalize
+
+def display_colorbar():
+    """Display a correct numeric colorbar for a shaded plot."""
+    y, x = np.mgrid[-4:2:200j, -4:2:200j]
+    z = 10 * np.cos(x**2 + y**2)
+
+    cmap = plt.cm.copper
+    ls = LightSource(315, 45)
+    rgb = ls.shade(z, cmap)
+
+    fig, ax = plt.subplots()
+    ax.imshow(rgb)
+
+    # Use a proxy artist for the colorbar...
+    im = ax.imshow(z, cmap=cmap)
+    im.remove()
+    fig.colorbar(im)
+
+    ax.set_title('Using a colorbar with a shaded plot', size='x-large')
+
+def avoid_outliers():
+    """Use a custom norm to control the displayed z-range of a shaded plot."""
+    y, x = np.mgrid[-4:2:200j, -4:2:200j]
+    z = 10 * np.cos(x**2 + y**2)
+
+    # Add some outliers...
+    z[100, 105] = 2000
+    z[120, 110] = -9000
+
+    ls = LightSource(315, 45)
+    fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(8, 4.5))
+
+    rgb = ls.shade(z, plt.cm.copper)
+    ax1.imshow(rgb)
+    ax1.set_title('Full range of data')
+
+    rgb = ls.shade(z, plt.cm.copper, vmin=-10, vmax=10)
+    ax2.imshow(rgb)
+    ax2.set_title('Manually set range')
+
+    fig.suptitle('Avoiding Outliers in Shaded Plots', size='x-large')
+
+def shade_other_data():
+    """Demonstrates displaying different variables through shade and color."""
+    y, x = np.mgrid[-4:2:200j, -4:2:200j]
+    z1 = np.sin(x**2) # Data to hillshade
+    z2 = np.cos(x**2 + y**2) # Data to color
+
+    norm=Normalize(z2.min(), z2.max())
+    cmap = plt.cm.jet
+
+    ls = LightSource(315, 45)
+    rgb = ls.shade_rgb(cmap(norm(z2)), z1)
+
+    fig, ax = plt.subplots()
+    ax.imshow(rgb)
+    ax.set_title('Shade by one variable, color by another', size='x-large')
+
+display_colorbar()
+avoid_outliers()
+shade_other_data()
+plt.show()
diff --git a/examples/specialty_plots/topographic_hillshading.py b/examples/specialty_plots/topographic_hillshading.py
@@ -0,0 +1,70 @@
+"""
+Demonstrates the visual effect of varying blend mode and vertical exaggeration
+on "hillshaded" plots.
+
+Note that the "overlay" and "soft" blend modes work well for complex surfaces
+such as this example, while the default "hsv" blend mode works best for smooth
+surfaces such as many mathematical functions.
+
+In most cases, hillshading is used purely for visual purposes, and *dx*/*dy*
+can be safely ignored. In that case, you can tweak *vert_exag* (vertical
+exaggeration) by trial and error to give the desired visual effect. However,
+this example demonstrates how to use the *dx* and *dy* kwargs to ensure that
+the *vert_exag* parameter is the true vertical exaggeration.
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.cbook import get_sample_data
+from matplotlib.colors import LightSource
+
+dem = np.load(get_sample_data('jacksboro_fault_dem.npz'))
+z = dem['elevation']
+
+#-- Optional dx and dy for accurate vertical exaggeration --------------------
+# If you need topographically accurate vertical exaggeration, or you don't want
+# to guess at what *vert_exag* should be, you'll need to specify the cellsize
+# of the grid (i.e. the *dx* and *dy* parameters).  Otherwise, any *vert_exag*
+# value you specify will be realitive to the grid spacing of your input data
+# (in other words, *dx* and *dy* default to 1.0, and *vert_exag* is calculated
+# relative to those parameters).  Similarly, *dx* and *dy* are assumed to be in
+# the same units as your input z-values.  Therefore, we'll need to convert the
+# given dx and dy from decimal degrees to meters.
+dx, dy = dem['dx'], dem['dy']
+dy = 111200 * dy
+dx = 111200 * dx * np.cos(np.radians(dem['ymin']))
+#-----------------------------------------------------------------------------
+
+# Shade from the northwest, with the sun 45 degrees from horizontal
+ls = LightSource(azdeg=315, altdeg=45)
+cmap = plt.cm.gist_earth
+
+fig, axes = plt.subplots(nrows=4, ncols=3, figsize=(8, 9))
+plt.setp(axes.flat, xticks=[], yticks=[])
+
+# Vary vertical exaggeration and blend mode and plot all combinations
+for col, ve in zip(axes.T, [0.1, 1, 10]):
+    # Show the hillshade intensity image in the first row
+    col[0].imshow(ls.hillshade(z, vert_exag=ve, dx=dx, dy=dy), cmap='gray')
+
+    # Place hillshaded plots with different blend modes in the rest of the rows
+    for ax, mode in zip(col[1:], ['hsv', 'overlay', 'soft']):
+        rgb = ls.shade(z, cmap=cmap, blend_mode=mode,
+                       vert_exag=ve, dx=dx, dy=dy)
+        ax.imshow(rgb)
+
+# Label rows and columns
+for ax, ve in zip(axes[0], [0.1, 1, 10]):
+    ax.set_title('{}'.format(ve), size=18)
+for ax, mode in zip(axes[:,0], ['Hillshade', 'hsv', 'overlay', 'soft']):
+    ax.set_ylabel(mode, size=18)
+
+# Group labels...
+axes[0,1].annotate('Vertical Exaggeration', (0.5, 1), xytext=(0, 30),
+                   textcoords='offset points', xycoords='axes fraction',
+                   ha='center', va='bottom', size=20)
+axes[2,0].annotate('Blend Mode', (0, 0.5), xytext=(-30, 0),
+                   textcoords='offset points', xycoords='axes fraction',
+                   ha='right', va='center', size=20, rotation=90)
+fig.subplots_adjust(bottom=0.05, right=0.95)
+
+plt.show()