diff --git a/doc/examples/plot_rag_draw.py b/doc/examples/plot_rag_draw.py
new file mode 100644
index 00000000..abf2ba08
--- /dev/null
+++ b/doc/examples/plot_rag_draw.py
@@ -0,0 +1,39 @@
+"""
+=====================================
+Drawing Region Adjacency Graphs (RAGs)
+======================================
+
+This example constructs a Region Adjacency Graph (RAG) and draws it with
+the `rag_draw` method.
+"""
+from skimage import graph, data, segmentation
+from matplotlib import pyplot as plt, colors
+
+
+img = data.coffee()
+labels = segmentation.slic(img, compactness=30, n_segments=400)
+g = graph.rag_mean_color(img, labels)
+
+out = graph.draw_rag(labels, g, img)
+plt.figure()
+plt.title("RAG with all edges shown in green.")
+plt.imshow(out)
+
+# The color palette used was taken from
+# http://www.colorcombos.com/color-schemes/2/ColorCombo2.html
+cmap = colors.ListedColormap(['#6599FF', '#ff9900'])
+out = graph.draw_rag(labels, g, img, node_color="#ffde00", colormap=cmap,
+                     thresh=30, desaturate=True)
+plt.figure()
+plt.title("RAG with edge weights less than 30, color "
+          "mapped between blue and orange.")
+plt.imshow(out)
+
+plt.figure()
+plt.title("All edges drawn with cubehelix colormap")
+cmap = plt.get_cmap('cubehelix')
+out = graph.draw_rag(labels, g, img, colormap=cmap,
+                     desaturate=True)
+
+plt.imshow(out)
+plt.show()
diff --git a/doc/release/release_dev.txt b/doc/release/release_dev.txt
index 1123eb77..cbeded27 100644
--- a/doc/release/release_dev.txt
+++ b/doc/release/release_dev.txt
@@ -14,9 +14,12 @@ http://scikit-image.org
 
 New Features
 ------------
-
-
-
+Region Adjacency Graphs
+ - Color distance RAGs (#1031)
+ - Threshold Cut on RAGs (#1031)
+ - Similarity RAGs (#1080)
+ - Normalized Cut on RAGs (#1080)
+ - RAG Drawing (#1087)
 
 Improvements
 ------------
diff --git a/skimage/graph/__init__.py b/skimage/graph/__init__.py
index aca73bc4..6da4fcc8 100644
--- a/skimage/graph/__init__.py
+++ b/skimage/graph/__init__.py
@@ -1,9 +1,10 @@
 from .spath import shortest_path
 from .mcp import MCP, MCP_Geometric, MCP_Connect, MCP_Flexible, route_through_array
-from .rag import rag_mean_color, RAG
 from .graph_cut import cut_threshold, cut_normalized
+from .rag import rag_mean_color, RAG, draw_rag
 ncut = cut_normalized
 
+
 __all__ = ['shortest_path',
            'MCP',
            'MCP_Geometric',
@@ -14,4 +15,5 @@ __all__ = ['shortest_path',
            'cut_threshold',
            'cut_normalized',
            'ncut',
+           'draw_rag',
            'RAG']
diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py
index d7710f0b..068e2150 100644
--- a/skimage/graph/graph_cut.py
+++ b/skimage/graph/graph_cut.py
@@ -128,7 +128,7 @@ def cut_normalized(labels, rag, thresh=0.001, num_cuts=10, in_place=True,
 
     _ncut_relabel(rag, thresh, num_cuts)
 
-    map_array = np.zeros(labels.max() + 1)
+    map_array = np.zeros(labels.max() + 1, dtype=labels.dtype)
     # Mapping from old labels to new
     for n, d in rag.nodes_iter(data=True):
         map_array[d['labels']] = d['ncut label']
diff --git a/skimage/graph/rag.py b/skimage/graph/rag.py
index df9b4ebf..5995fb5d 100644
--- a/skimage/graph/rag.py
+++ b/skimage/graph/rag.py
@@ -9,10 +9,17 @@ except ImportError:
                 raise ImportError(msg)
     import warnings
     warnings.warn(msg)
+
 import numpy as np
 from scipy.ndimage import filters
 from scipy import ndimage as nd
 import math
+from .. import draw, measure, segmentation, util, color
+try:
+    from matplotlib import colors
+    from matplotlib import cm
+except ImportError:
+    pass
 
 
 def min_weight(graph, src, dst, n):
@@ -237,3 +244,110 @@ def rag_mean_color(image, labels, connectivity=2, mode='distance',
             raise ValueError("The mode '%s' is not recognised" % mode)
 
     return graph
+
+
+def draw_rag(labels, rag, img, border_color=None, node_color='#ffff00',
+             edge_color='#00ff00', colormap=None, thresh=np.inf,
+             desaturate=False, in_place=True):
+    """Draw a Region Adjacency Graph on an image.
+
+    Given a labelled image and its corresponding RAG, draw the nodes and edges
+    of the RAG on the image with the specified colors. Nodes are marked by
+    the centroids of the corresponding regions.
+
+    Parameters
+    ----------
+    labels : ndarray, shape (M, N)
+        The labelled image.
+    rag : RAG
+        The Region Adjacency Graph.
+    img : ndarray, shape (M, N, 3)
+        Input image.
+    border_color : colorspec, optional
+        Any matplotlib colorspec.
+    node_color : colorspec, optional
+        Any matplotlib colorspec. Yellow by default.
+    edge_color : colorspec, optional
+        Any matplotlib colorspec. Green by default.
+    colormap : colormap, optional
+        Any matplotlib colormap. If specified the edges are colormapped with
+        the specified color map.
+    thresh : float, optional
+        Edges with weight below `thresh` are not drawn, or considered for color
+        mapping.
+    desaturate : bool, optional
+        Convert the image to grayscale before displaying. Particularly helps
+        visualization when using the `colormap` option.
+    in_place : bool, optional
+        If set, the RAG is modified in place. For each node `n` the function
+        will set a new attribute ``rag.node[n]['centroid']``.
+
+    Returns
+    -------
+    out : ndarray, shape (M, N, 3)
+        The image with the RAG drawn.
+
+    Examples
+    --------
+    >>> from skimage import data, graph, segmentation
+    >>> img = data.coffee()
+    >>> labels = segmentation.slic(img)
+    >>> g =  graph.rag_mean_color(img, labels)
+    >>> out = graph.draw_rag(labels, g, img)
+    """
+    if not in_place:
+        rag = rag.copy()
+
+    if desaturate:
+        img = color.rgb2gray(img)
+        img = color.gray2rgb(img)
+
+    out = util.img_as_float(img, force_copy=True)
+    cc = colors.ColorConverter()
+
+    edge_color = cc.to_rgb(edge_color)
+    node_color = cc.to_rgb(node_color)
+
+    # Handling the case where one node has multiple labels
+    # offset is 1 so that regionprops does not ignore 0
+    offset = 1
+    map_array = np.arange(labels.max() + 1)
+    for n, d in rag.nodes_iter(data=True):
+        for label in d['labels']:
+            map_array[label] = offset
+        offset += 1
+
+    rag_labels = map_array[labels]
+    regions = measure.regionprops(rag_labels)
+
+    for (n, data), region in zip(rag.nodes_iter(data=True), regions):
+        data['centroid'] = region['centroid']
+
+    if border_color is not None:
+        border_color = cc.to_rgb(border_color)
+        out = segmentation.mark_boundaries(out, rag_labels, color=border_color)
+
+    if colormap is not None:
+        edge_weight_list = [d['weight'] for x, y, d in
+                            rag.edges_iter(data=True) if d['weight'] < thresh]
+        norm = colors.Normalize()
+        norm.autoscale(edge_weight_list)
+        smap = cm.ScalarMappable(norm, colormap)
+
+    for n1, n2, data in rag.edges_iter(data=True):
+
+        if data['weight'] >= thresh:
+            continue
+        r1, c1 = map(int, rag.node[n1]['centroid'])
+        r2, c2 = map(int, rag.node[n2]['centroid'])
+        line = draw.line(r1, c1, r2, c2)
+
+        if colormap is not None:
+            out[line] = smap.to_rgba([data['weight']])[0][:-1]
+        else:
+            out[line] = edge_color
+
+        circle = draw.circle(r1, c1, 2)
+        out[circle] = node_color
+
+    return out