diff --git a/TODO.txt b/TODO.txt
index a96b3731..1af04a20 100644
--- a/TODO.txt
+++ b/TODO.txt
@@ -16,6 +16,10 @@ Version 0.13
 * Remove deprecated edge filters `hsobel`, `vsobel`, `hscharr`, `vscharr`,
   `hprewitt`, `vprewitt`, `roberts_positive_diagonal`,
   `roberts_negative_diagonal` in `skimage/filters/edges.py`
+* Remove supported for renamed edge mode, 'nearest' (it is now 'edge').  This
+  involves removing the function _mode_deprecations from skimage._shared.utils
+  as well as any uses of _mode_deprecations from restoration/_denoise.py,
+  _shared/interpolation.pyx, transform/_geometric.py, and transform/_warps.py
 
 
 
diff --git a/bento.info b/bento.info
index 659247d1..1d7f3a28 100644
--- a/bento.info
+++ b/bento.info
@@ -107,6 +107,9 @@ Library:
     Extension: skimage._shared.transform
         Sources:
             skimage/_shared/transform.pyx
+    Extension: skimage._shared.interpolation
+        Sources:
+            skimage/_shared/interpolation.pyx
     Extension: skimage.segmentation._slic
         Sources:
             skimage/segmentation/_slic.pyx
diff --git a/doc/examples/plot_edge_modes.py b/doc/examples/plot_edge_modes.py
new file mode 100644
index 00000000..d68abf02
--- /dev/null
+++ b/doc/examples/plot_edge_modes.py
@@ -0,0 +1,37 @@
+"""
+=========================
+Interpolation: Edge Modes
+=========================
+
+This example illustrates the different edge modes available during
+interpolation in routines such as ``skimage.transform.rescale`` and
+``skimage.transform.resize``.
+"""
+from skimage._shared.interpolation import extend_image
+import skimage.data
+import matplotlib.pyplot as plt
+import numpy as np
+
+img = np.zeros((16, 16))
+img[:8, :8] += 1
+img[:4, :4] += 1
+img[:2, :2] += 1
+img[:1, :1] += 2
+img[8, 8] = 4
+
+modes = ['constant', 'edge', 'wrap', 'reflect', 'symmetric']
+fig, axes = plt.subplots(1, 5, figsize=(15, 5))
+for n, mode in enumerate(modes):
+    img_extended = extend_image(img, pad=img.shape[0], mode=mode)
+    axes[n].imshow(img_extended, cmap=plt.cm.gray, interpolation='nearest')
+    axes[n].plot([15.5, 15.5], [15.5, 31.5], 'y--', linewidth=0.5)
+    axes[n].plot([31.5, 31.5], [15.5, 31.5], 'y--', linewidth=0.5)
+    axes[n].plot([15.5, 31.5], [15.5, 15.5], 'y--', linewidth=0.5)
+    axes[n].plot([15.5, 31.5], [31.5, 31.5], 'y--', linewidth=0.5)
+    axes[n].set_axis_off()
+    axes[n].set_aspect('equal')
+    axes[n].set_title(mode)
+
+plt.tight_layout()
+
+plt.show()
diff --git a/skimage/_shared/interpolation.pxd b/skimage/_shared/interpolation.pxd
index 3b4ac538..aa58684d 100644
--- a/skimage/_shared/interpolation.pxd
+++ b/skimage/_shared/interpolation.pxd
@@ -2,6 +2,20 @@
 #cython: boundscheck=False
 #cython: nonecheck=False
 #cython: wraparound=False
+"""
+Note: All edge modes implemented here follow the corresponding numpy.pad
+conventions.
+
+The table below illustrates the behavior for the array [1, 2, 3, 4], if padded
+by 4 values on each side:
+
+                               pad     original    pad
+    constant (with c=0) :    0 0 0 0 | 1 2 3 4 | 0 0 0 0
+    wrap                :    1 2 3 4 | 1 2 3 4 | 1 2 3 4
+    symmetric           :    4 3 2 1 | 1 2 3 4 | 4 3 2 1
+    edge                :    1 1 1 1 | 1 2 3 4 | 4 4 4 4
+    reflect             :    3 4 3 2 | 1 2 3 4 | 3 2 1 2
+"""
 from libc.math cimport ceil, floor
 
 
@@ -24,8 +38,8 @@ cdef inline double nearest_neighbour_interpolation(double* image,
         Shape of image.
     r, c : double
         Position at which to interpolate.
-    mode : {'C', 'W', 'R', 'N'}
-        Wrapping mode. Constant, Wrap, Reflect or Nearest.
+    mode : {'C', 'W', 'S', 'E', 'R'}
+        Wrapping mode. Constant, Wrap, Symmetric, Edge or Reflect.
     cval : double
         Constant value to use for constant mode.
 
@@ -52,8 +66,8 @@ cdef inline double bilinear_interpolation(double* image, Py_ssize_t rows,
         Shape of image.
     r, c : double
         Position at which to interpolate.
-    mode : {'C', 'W', 'R', 'N'}
-        Wrapping mode. Constant, Wrap, Reflect or Nearest.
+    mode : {'C', 'W', 'S', 'E', 'R'}
+        Wrapping mode. Constant, Wrap, Symmetric, Edge or Reflect.
     cval : double
         Constant value to use for constant mode.
 
@@ -119,8 +133,8 @@ cdef inline double biquadratic_interpolation(double* image, Py_ssize_t rows,
         Shape of image.
     r, c : double
         Position at which to interpolate.
-    mode : {'C', 'W', 'R', 'N'}
-        Wrapping mode. Constant, Wrap, Reflect or Nearest.
+    mode : {'C', 'W', 'S', 'E', 'R'}
+        Wrapping mode. Constant, Wrap, Symmetric, Edge or Reflect.
     cval : double
         Constant value to use for constant mode.
 
@@ -192,8 +206,8 @@ cdef inline double bicubic_interpolation(double* image, Py_ssize_t rows,
         Shape of image.
     r, c : double
         Position at which to interpolate.
-    mode : {'C', 'W', 'R', 'N'}
-        Wrapping mode. Constant, Wrap, Reflect or Nearest.
+    mode : {'C', 'W', 'S', 'E', 'R'}
+        Wrapping mode. Constant, Wrap, Symmetric, Edge or Reflect.
     cval : double
         Constant value to use for constant mode.
 
@@ -248,8 +262,8 @@ cdef inline double get_pixel2d(double* image, Py_ssize_t rows, Py_ssize_t cols,
         Shape of image.
     r, c : int
         Position at which to get the pixel.
-    mode : {'C', 'W', 'R', 'N'}
-        Wrapping mode. Constant, Wrap, Reflect or Nearest.
+    mode : {'C', 'W', 'S', 'E', 'R'}
+        Wrapping mode. Constant, Wrap, Symmetric, Edge or Reflect.
     cval : double
         Constant value to use for constant mode.
 
@@ -281,8 +295,8 @@ cdef inline double get_pixel3d(double* image, Py_ssize_t rows, Py_ssize_t cols,
         Shape of image.
     r, c, d : int
         Position at which to get the pixel.
-    mode : {'C', 'W', 'R', 'N'}
-        Wrapping mode. Constant, Wrap, Reflect or Nearest.
+    mode : {'C', 'W', 'S', 'E', 'R'}
+        Wrapping mode. Constant, Wrap, Symmetric, Edge or Reflect.
     cval : double
         Constant value to use for constant mode.
 
@@ -312,33 +326,40 @@ cdef inline Py_ssize_t coord_map(Py_ssize_t dim, long coord, char mode) nogil:
         Maximum coordinate.
     coord : int
         Coord provided by user.  May be < 0 or > dim.
-    mode : {'W', 'R', 'N'}
-        Whether to wrap or reflect the coordinate if it
-        falls outside [0, dim).
-
+    mode : {'W', 'S', 'R', 'E'}
+        Whether to wrap, symmetric reflect, reflect or use the nearest
+        coordinate if `coord` falls outside [0, dim).
     """
-    dim = dim - 1
-    if mode == 'R': # reflect
+    cdef Py_ssize_t cmax
+    cmax = dim - 1
+    if mode == 'S': # symmetric
         if coord < 0:
-            # How many times times does the coordinate wrap?
-            if <Py_ssize_t>(-coord / dim) % 2 != 0:
-                return dim - <Py_ssize_t>(-coord % dim)
-            else:
-                return <Py_ssize_t>(-coord % dim)
-        elif coord > dim:
+            coord = -coord - 1
+        if coord > cmax:
             if <Py_ssize_t>(coord / dim) % 2 != 0:
-                return <Py_ssize_t>(dim - (coord % dim))
+                return <Py_ssize_t>(cmax - (coord % dim))
             else:
                 return <Py_ssize_t>(coord % dim)
     elif mode == 'W': # wrap
         if coord < 0:
-            return <Py_ssize_t>(dim - (-coord % dim))
-        elif coord > dim:
+            return <Py_ssize_t>(cmax - ((-coord - 1) % dim))
+        elif coord > cmax:
             return <Py_ssize_t>(coord % dim)
-    elif mode == 'N': # nearest
+    elif mode == 'E': # edge
         if coord < 0:
             return 0
-        elif coord > dim:
-            return dim
-
+        elif coord > cmax:
+            return cmax
+    elif mode == 'R': # reflect (mirror)
+        if coord < 0:
+            # How many times times does the coordinate wrap?
+            if <Py_ssize_t>(-coord / cmax) % 2 != 0:
+                return cmax - <Py_ssize_t>(-coord % cmax)
+            else:
+                return <Py_ssize_t>(-coord % cmax)
+        elif coord > cmax:
+            if <Py_ssize_t>(coord / cmax) % 2 != 0:
+                return <Py_ssize_t>(cmax - (coord % cmax))
+            else:
+                return <Py_ssize_t>(coord % cmax)
     return coord
diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx
new file mode 100644
index 00000000..583fe43c
--- /dev/null
+++ b/skimage/_shared/interpolation.pyx
@@ -0,0 +1,57 @@
+from interpolation cimport coord_map, get_pixel2d
+import numpy as np
+cimport numpy as cnp
+from .utils import _mode_deprecations
+
+
+def coord_map_py(Py_ssize_t dim, long coord, mode):
+    """interpolation.coord_map python wrapper"""
+    cdef char mode_c = ord(mode[0].upper())
+    return coord_map(dim, coord, mode_c)
+
+
+def extend_image(image, pad=10, mode='constant', cval=0):
+    """Pad a 2D image by ``pad`` pixels on each side.
+
+    Parameters
+    ----------
+    image : ndarray
+        Input image.
+    pad : int, optional
+        The number of pixels to pad around the border
+    mode : {'constant', 'edge', 'symmetric', 'reflect', 'wrap'}, optional
+        Points outside the boundaries of the input are filled according
+        to the given mode.
+    cval : float, optional
+        Used in conjunction with mode 'constant', the value outside
+        the image boundaries.
+
+    Returns
+    -------
+    extended : ndarray
+        The extended version of the input image.
+
+    Note
+    ----
+    For image padding, ``skimage.util.pad`` should be used instead.  This
+    function is intended only for testing get_pixel2d and demonstrating the
+    coordinate mapping modes implemented in ``coord_map``.
+    """
+    mode = _mode_deprecations(mode)
+    cdef:
+        Py_ssize_t rows = image.shape[0]
+        Py_ssize_t cols = image.shape[1]
+        long ro, co
+        char mode_c = ord(mode[0].upper())
+
+    image = np.ascontiguousarray(image.astype(np.float64))
+    output_shape = np.asarray(image.shape) + 2 * pad
+    extended = np.zeros(output_shape, dtype=image.dtype)
+    for r in range(-pad, rows + pad):
+        for c in range(-pad, cols + pad):
+            ro = r + pad
+            co = c + pad
+            extended[ro, co] = get_pixel2d(<double*> cnp.PyArray_DATA(image),
+                                           rows, cols, <long> r, <long> c,
+                                           mode_c, <double> cval)
+    return extended
diff --git a/skimage/_shared/setup.py b/skimage/_shared/setup.py
index 4c36b97e..066a856f 100644
--- a/skimage/_shared/setup.py
+++ b/skimage/_shared/setup.py
@@ -15,11 +15,12 @@ def configuration(parent_package='', top_path=None):
 
     cython(['geometry.pyx'], working_path=base_path)
     cython(['transform.pyx'], working_path=base_path)
+    cython(['interpolation.pyx'], working_path=base_path)
 
     config.add_extension('geometry', sources=['geometry.c'])
     config.add_extension('transform', sources=['transform.c'],
                          include_dirs=[get_numpy_include_dirs()])
-
+    config.add_extension('interpolation', sources=['interpolation.c'])
     return config
 
 
diff --git a/skimage/_shared/tests/test_interpolation.py b/skimage/_shared/tests/test_interpolation.py
new file mode 100644
index 00000000..61a24ea9
--- /dev/null
+++ b/skimage/_shared/tests/test_interpolation.py
@@ -0,0 +1,27 @@
+from skimage._shared.interpolation import coord_map_py
+from numpy.testing import assert_array_equal
+
+
+def test_coord_map():
+    symmetric = [coord_map_py(4, n, 'S') for n in range(-6, 6)]
+    expected_symmetric = [2, 3, 3, 2, 1, 0, 0, 1, 2, 3, 3, 2]
+    assert_array_equal(symmetric, expected_symmetric)
+
+    wrap = [coord_map_py(4, n, 'W') for n in range(-6, 6)]
+    expected_wrap = [2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1]
+    assert_array_equal(wrap, expected_wrap)
+
+    edge = [coord_map_py(4, n, 'E') for n in range(-6, 6)]
+    expected_edge = [0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 3, 3]
+    assert_array_equal(edge, expected_edge)
+
+    reflect = [coord_map_py(4, n, 'R') for n in range(-6, 6)]
+    expected_reflect = [0, 1, 2, 3, 2, 1, 0, 1, 2, 3, 2, 1]
+    assert_array_equal(reflect, expected_reflect)
+
+    constant = [coord_map_py(4, n, 'C') for n in range(-6, 6)]
+    expected_constant = [0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 0, 0]
+    assert_array_equal(constant, expected_constant)
+
+    other = [coord_map_py(4, n, 'undefined') for n in range(-6, 6)]
+    assert_array_equal(other, list(range(-6, 6)))
diff --git a/skimage/_shared/utils.py b/skimage/_shared/utils.py
index 43fda35d..eadc3015 100644
--- a/skimage/_shared/utils.py
+++ b/skimage/_shared/utils.py
@@ -163,3 +163,14 @@ def assert_nD(array, ndim, arg_name='image'):
         ndim = [ndim]
     if not array.ndim in ndim:
         raise ValueError(msg % (arg_name, '-or-'.join([str(n) for n in ndim])))
+
+
+def _mode_deprecations(mode):
+    """ to be used by functions to update deprecated mode names in
+    `skimage._shared.interpolation.pyx`."""
+    if mode.lower() == 'nearest':
+        warnings.warn(skimage_deprecation(
+            "Mode 'nearest' has been renamed 'edge'. Mode 'nearest' will be "
+            "removed in a future release."))
+        mode = 'edge'
+    return mode
diff --git a/skimage/filters/_gabor.py b/skimage/filters/_gabor.py
index 3eeccc0f..4d249d2c 100644
--- a/skimage/filters/_gabor.py
+++ b/skimage/filters/_gabor.py
@@ -129,7 +129,7 @@ def gabor_filter(image, frequency, theta=0, bandwidth=1, sigma_x=None,
         deviations.
     offset : float, optional
         Phase offset of harmonic function in radians.
-    mode : string, optional
+    mode : {'constant', 'nearest', 'reflect', 'mirror', 'wrap'}, optional
         Mode used to convolve image with a kernel, passed to `ndi.convolve`
     cval : scalar, optional
         Value to fill past edges of input if `mode` of convolution is
diff --git a/skimage/measure/profile.py b/skimage/measure/profile.py
index 819b8cff..9dbc2d7a 100644
--- a/skimage/measure/profile.py
+++ b/skimage/measure/profile.py
@@ -21,7 +21,7 @@ def profile_line(img, src, dst, linewidth=1,
     order : int in {0, 1, 2, 3, 4, 5}, optional
         The order of the spline interpolation to compute image values at
         non-integer coordinates. 0 means nearest-neighbor interpolation.
-    mode : string, one of {'constant', 'nearest', 'reflect', 'wrap'}, optional
+    mode : {'constant', 'nearest', 'reflect', 'mirror', 'wrap'}, optional
         How to compute any values falling outside of the image.
     cval : float, optional
         If `mode` is 'constant', what constant value to use outside the image.
diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 10de6c68..c44ab1f1 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -2,6 +2,7 @@
 import numpy as np
 from .. import img_as_float
 from ..restoration._denoise_cy import _denoise_bilateral, _denoise_tv_bregman
+from .._shared.utils import _mode_deprecations
 
 
 def denoise_bilateral(image, win_size=5, sigma_range=None, sigma_spatial=1,
@@ -37,9 +38,9 @@ def denoise_bilateral(image, win_size=5, sigma_range=None, sigma_spatial=1,
     bins : int
         Number of discrete values for gaussian weights of color filtering.
         A larger value results in improved accuracy.
-    mode : string
+    mode : {'constant', 'edge', 'symmetric', 'reflect', 'wrap'}
         How to handle values outside the image borders. See
-        `scipy.ndimage.map_coordinates` for detail.
+        `numpy.pad` for detail.
     cval : string
         Used in conjunction with mode 'constant', the value outside
         the image boundaries.
@@ -54,6 +55,7 @@ def denoise_bilateral(image, win_size=5, sigma_range=None, sigma_spatial=1,
     .. [1] http://users.soe.ucsc.edu/~manduchi/Papers/ICCV98.pdf
 
     """
+    mode = _mode_deprecations(mode)
     return _denoise_bilateral(image, win_size, sigma_range, sigma_spatial,
                               bins, mode, cval)
 
diff --git a/skimage/restoration/_denoise_cy.pyx b/skimage/restoration/_denoise_cy.pyx
index 7d3a82b9..b679b488 100644
--- a/skimage/restoration/_denoise_cy.pyx
+++ b/skimage/restoration/_denoise_cy.pyx
@@ -105,9 +105,9 @@ def _denoise_bilateral(image, Py_ssize_t win_size, sigma_range,
     centres = <double*>malloc(dims * sizeof(double))
     total_values = <double*>malloc(dims * sizeof(double))
 
-    if mode not in ('constant', 'wrap', 'reflect', 'nearest'):
-        raise ValueError("Invalid mode specified.  Please use "
-                         "`constant`, `nearest`, `wrap` or `reflect`.")
+    if mode not in ('constant', 'wrap', 'symmetric', 'reflect', 'edge'):
+        raise ValueError("Invalid mode specified.  Please use `constant`, "
+                         "`edge`, `wrap`, `symmetric` or `reflect`.")
     cdef char cmode = ord(mode[0].upper())
 
     for r in range(rows):
diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 1245a127..3d459b36 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -5,11 +5,24 @@ import numpy as np
 from scipy import spatial
 from scipy import ndimage as ndi
 
-from .._shared.utils import get_bound_method_class, safe_as_int
+from .._shared.utils import (get_bound_method_class, safe_as_int,
+                             _mode_deprecations)
 from ..util import img_as_float
+
 from ._warps_cy import _warp_fast
 
 
+def _to_ndimage_mode(mode):
+    """ Convert from a numpy.pad mode name to the corresponding ndimage
+    mode. """
+    mode = _mode_deprecations(mode.lower())
+    mode_translation_dict = dict(edge='nearest', symmetric='reflect',
+                                 reflect='mirror')
+    if mode in mode_translation_dict:
+        mode = mode_translation_dict[mode]
+    return mode
+
+
 def _center_and_normalize_points(points):
     """Center and normalize image points.
 
@@ -1128,9 +1141,9 @@ def _clip_warp_output(input_image, output_image, order, mode, cval, clip):
     order : int, optional
         The order of the spline interpolation, default is 1. The order has to
         be in the range 0-5. See `skimage.transform.warp` for detail.
-    mode : string, optional
+    mode : {'constant', 'edge', 'symmetric', 'reflect', 'wrap'}, optional
         Points outside the boundaries of the input are filled according
-        to the given mode ('constant', 'nearest', 'reflect' or 'wrap').
+        to the given mode.  Modes match the behaviour of `numpy.pad`.
     cval : float, optional
         Used in conjunction with mode 'constant', the value outside
         the image boundaries.
@@ -1140,7 +1153,6 @@ def _clip_warp_output(input_image, output_image, order, mode, cval, clip):
         produce values outside the given input range.
 
     """
-
     if clip and order != 0:
         min_val = input_image.min()
         max_val = input_image.max()
@@ -1211,9 +1223,9 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
          - 3: Bi-cubic
          - 4: Bi-quartic
          - 5: Bi-quintic
-    mode : string, optional
+    mode : {'constant', 'edge', 'symmetric', 'reflect', 'wrap'}, optional
         Points outside the boundaries of the input are filled according
-        to the given mode ('constant', 'nearest', 'reflect' or 'wrap').
+        to the given mode.  Modes match the behaviour of `numpy.pad`.
     cval : float, optional
         Used in conjunction with mode 'constant', the value outside
         the image boundaries.
@@ -1294,7 +1306,7 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
     >>> warped = warp(cube, coords)
 
     """
-
+    mode = _mode_deprecations(mode)
     image = _convert_warp_input(image, preserve_range)
 
     input_shape = np.array(image.shape)
@@ -1388,8 +1400,9 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
         # Pre-filtering not necessary for order 0, 1 interpolation
         prefilter = order > 1
 
+        ndi_mode = _to_ndimage_mode(mode)
         warped = ndi.map_coordinates(image, coords, prefilter=prefilter,
-                                     mode=mode, order=order, cval=cval)
+                                     mode=ndi_mode, order=order, cval=cval)
 
 
     _clip_warp_output(image, warped, order, mode, cval, clip)
diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py
index f069122a..2f2cf553 100644
--- a/skimage/transform/_warps.py
+++ b/skimage/transform/_warps.py
@@ -3,7 +3,8 @@ from scipy import ndimage as ndi
 
 from ..measure import block_reduce
 from ._geometric import (warp, SimilarityTransform, AffineTransform,
-                         _convert_warp_input, _clip_warp_output)
+                         _convert_warp_input, _clip_warp_output,
+                         _to_ndimage_mode)
 
 
 def resize(image, output_shape, order=1, mode='constant', cval=0, clip=True,
@@ -35,9 +36,9 @@ def resize(image, output_shape, order=1, mode='constant', cval=0, clip=True,
     order : int, optional
         The order of the spline interpolation, default is 1. The order has to
         be in the range 0-5. See `skimage.transform.warp` for detail.
-    mode : string, optional
+    mode : {'constant', 'edge', 'symmetric', 'reflect', 'wrap'}, optional
         Points outside the boundaries of the input are filled according
-        to the given mode ('constant', 'nearest', 'reflect' or 'wrap').
+        to the given mode.  Modes match the behaviour of `numpy.pad`.
     cval : float, optional
         Used in conjunction with mode 'constant', the value outside
         the image boundaries.
@@ -49,6 +50,14 @@ def resize(image, output_shape, order=1, mode='constant', cval=0, clip=True,
         Whether to keep the original range of values. Otherwise, the input
         image is converted according to the conventions of `img_as_float`.
 
+    Note
+    ----
+    Modes 'reflect' and 'symmetric' are similar, but differ in whether the edge
+    pixels are duplicated during the reflection.  As an example, if an array
+    has values [0, 1, 2] and was padded to the right by four values using
+    symmetric, the result would be [0, 1, 2, 2, 1, 0, 0], while for reflect it
+    would be [0, 1, 2, 1, 0, 1, 2].
+
     Examples
     --------
     >>> from skimage import data
@@ -68,8 +77,11 @@ def resize(image, output_shape, order=1, mode='constant', cval=0, clip=True,
     # 3-dimensional interpolation
     if len(output_shape) == 3 and (image.ndim == 2
                                    or output_shape[2] != image.shape[2]):
+        ndi_mode = _to_ndimage_mode(mode)
         dim = output_shape[2]
-        orig_dim = 1 if image.ndim == 2 else image.shape[2]
+        if image.ndim == 2:
+            image = image[:, :, np.newaxis]
+        orig_dim = image.shape[2]
         dim_scale = float(orig_dim) / dim
 
         map_rows, map_cols, map_dims = np.mgrid[:rows, :cols, :dim]
@@ -82,7 +94,7 @@ def resize(image, output_shape, order=1, mode='constant', cval=0, clip=True,
         image = _convert_warp_input(image, preserve_range)
 
         out = ndi.map_coordinates(image, coord_map, order=order,
-                                  mode=mode, cval=cval)
+                                  mode=ndi_mode, cval=cval)
 
         _clip_warp_output(image, out, order, mode, cval, clip)
 
@@ -136,9 +148,9 @@ def rescale(image, scale, order=1, mode='constant', cval=0, clip=True,
     order : int, optional
         The order of the spline interpolation, default is 1. The order has to
         be in the range 0-5. See `skimage.transform.warp` for detail.
-    mode : string, optional
+    mode : {'constant', 'edge', 'symmetric', 'reflect', 'wrap'}, optional
         Points outside the boundaries of the input are filled according
-        to the given mode ('constant', 'nearest', 'reflect' or 'wrap').
+        to the given mode.  Modes match the behaviour of `numpy.pad`.
     cval : float, optional
         Used in conjunction with mode 'constant', the value outside
         the image boundaries.
@@ -204,9 +216,9 @@ def rotate(image, angle, resize=False, center=None, order=1, mode='constant',
     order : int, optional
         The order of the spline interpolation, default is 1. The order has to
         be in the range 0-5. See `skimage.transform.warp` for detail.
-    mode : string, optional
+    mode : {'constant', 'edge', 'symmetric', 'reflect', 'wrap'}, optional
         Points outside the boundaries of the input are filled according
-        to the given mode ('constant', 'nearest', 'reflect' or 'wrap').
+        to the given mode.  Modes match the behaviour of `numpy.pad`.
     cval : float, optional
         Used in conjunction with mode 'constant', the value outside
         the image boundaries.
@@ -358,9 +370,9 @@ def swirl(image, center=None, strength=1, radius=100, rotation=0,
     order : int, optional
         The order of the spline interpolation, default is 1. The order has to
         be in the range 0-5. See `skimage.transform.warp` for detail.
-    mode : string, optional
+    mode : {'constant', 'edge', 'symmetric', 'reflect', 'wrap'}, optional
         Points outside the boundaries of the input are filled according
-        to the given mode ('constant', 'nearest', 'reflect' or 'wrap').
+        to the given mode.  Modes match the behaviour of `numpy.pad`.
     cval : float, optional
         Used in conjunction with mode 'constant', the value outside
         the image boundaries.
diff --git a/skimage/transform/_warps_cy.pyx b/skimage/transform/_warps_cy.pyx
index aa2f3b6e..a3167016 100644
--- a/skimage/transform/_warps_cy.pyx
+++ b/skimage/transform/_warps_cy.pyx
@@ -70,20 +70,29 @@ def _warp_fast(cnp.ndarray image, cnp.ndarray H, output_shape=None,
         * 1: Bi-linear (default)
         * 2: Bi-quadratic
         * 3: Bi-cubic
-    mode : {'constant', 'reflect', 'wrap', 'nearest'}, optional
-        How to handle values outside the image borders (default is constant).
+    mode : {'constant', 'edge', 'symmetric', 'reflect', 'wrap'}, optional
+        Points outside the boundaries of the input are filled according
+        to the given mode.  Modes match the behaviour of `numpy.pad`.
     cval : string, optional (default 0)
         Used in conjunction with mode 'C' (constant), the value
         outside the image boundaries.
 
+    Note
+    ----
+    Modes 'reflect' and 'symmetric' are similar, but differ in whether the edge
+    pixels are duplicated during the reflection.  As an example, if an array
+    has values [0, 1, 2] and was padded to the right by four values using
+    symmetric, the result would be [0, 1, 2, 2, 1, 0, 0], while for reflect it
+    would be [0, 1, 2, 1, 0, 1, 2].
+
     """
 
     cdef double[:, ::1] img = np.ascontiguousarray(image, dtype=np.double)
     cdef double[:, ::1] M = np.ascontiguousarray(H)
 
-    if mode not in ('constant', 'wrap', 'reflect', 'nearest'):
-        raise ValueError("Invalid mode specified.  Please use "
-                         "`constant`, `nearest`, `wrap` or `reflect`.")
+    if mode not in ('constant', 'wrap', 'symmetric', 'reflect', 'edge'):
+        raise ValueError("Invalid mode specified.  Please use `constant`, "
+                         "`edge`, `wrap`, `reflect` or `symmetric`.")
     cdef char mode_c = ord(mode[0].upper())
 
     cdef Py_ssize_t out_r, out_c
diff --git a/skimage/transform/pyramids.py b/skimage/transform/pyramids.py
index fe2d26d1..958e9d65 100644
--- a/skimage/transform/pyramids.py
+++ b/skimage/transform/pyramids.py
@@ -45,7 +45,7 @@ def pyramid_reduce(image, downscale=2, sigma=None, order=1,
     order : int, optional
         Order of splines used in interpolation of downsampling. See
         `skimage.transform.warp` for detail.
-    mode : {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}, optional
+    mode : {'reflect', 'constant', 'edge', 'symmetric', 'wrap'}, optional
         The mode parameter determines how the array borders are handled, where
         cval is the value when mode is equal to 'constant'.
     cval : float, optional
@@ -99,7 +99,7 @@ def pyramid_expand(image, upscale=2, sigma=None, order=1,
     order : int, optional
         Order of splines used in interpolation of upsampling. See
         `skimage.transform.warp` for detail.
-    mode : {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}, optional
+    mode : {'reflect', 'constant', 'edge', 'symmetric', 'wrap'}, optional
         The mode parameter determines how the array borders are handled, where
         cval is the value when mode is equal to 'constant'.
     cval : float, optional
@@ -164,7 +164,7 @@ def pyramid_gaussian(image, max_layer=-1, downscale=2, sigma=None, order=1,
     order : int, optional
         Order of splines used in interpolation of downsampling. See
         `skimage.transform.warp` for detail.
-    mode : {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}, optional
+    mode : {'reflect', 'constant', 'edge', 'symmetric', 'wrap'}, optional
         The mode parameter determines how the array borders are handled, where
         cval is the value when mode is equal to 'constant'.
     cval : float, optional
@@ -245,7 +245,7 @@ def pyramid_laplacian(image, max_layer=-1, downscale=2, sigma=None, order=1,
     order : int, optional
         Order of splines used in interpolation of downsampling. See
         `skimage.transform.warp` for detail.
-    mode : {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}, optional
+    mode : {'reflect', 'constant', 'edge', 'symmetric', 'wrap'}, optional
         The mode parameter determines how the array borders are handled, where
         cval is the value when mode is equal to 'constant'.
     cval : float, optional
diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py
index 4de59359..a8dfcb87 100644
--- a/skimage/transform/tests/test_warps.py
+++ b/skimage/transform/tests/test_warps.py
@@ -181,6 +181,16 @@ def test_resize3d_resize():
     assert_almost_equal(resized, ref)
 
 
+def test_resize3d_2din_3dout():
+    # 3D output with 2D input
+    x = np.zeros((5, 5), dtype=np.double)
+    x[1, 1] = 1
+    resized = resize(x, (10, 10, 1), order=0)
+    ref = np.zeros((10, 10, 1))
+    ref[2:4, 2:4] = 1
+    assert_almost_equal(resized, ref)
+
+
 def test_resize3d_bilinear():
     # bilinear 3rd dimension
     x = np.zeros((5, 5, 2), dtype=np.double)
diff --git a/skimage/util/apply_parallel.py b/skimage/util/apply_parallel.py
index fd50ba56..c19b4bfd 100644
--- a/skimage/util/apply_parallel.py
+++ b/skimage/util/apply_parallel.py
@@ -68,13 +68,18 @@ def apply_parallel(function, array, chunks=None, depth=0, mode=None,
     depth : int, optional
         Integer equal to the depth of the added boundary cells. Defaults to
         zero.
-    mode : 'reflect', 'periodic', 'wrap', 'nearest', optional
-        type of external boundary padding
+    mode : {'reflect', 'symmetric', 'periodic', 'wrap', 'nearest', 'edge'}, optional
+        type of external boundary padding.
     extra_arguments : tuple, optional
         Tuple of arguments to be passed to the function.
     extra_keywords : dictionary, optional
         Dictionary of keyword arguments to be passed to the function.
 
+    Notes
+    -----
+    Numpy edge modes `symmetric`, `wrap` and `edge` are converted to the
+    equivalent `dask` boundary modes `reflect`, `periodic` and `nearest`,
+    respectively.
     """
     import dask.array as da
 
@@ -88,6 +93,10 @@ def apply_parallel(function, array, chunks=None, depth=0, mode=None,
 
     if mode == 'wrap':
         mode = 'periodic'
+    elif mode == 'symmetric':
+        mode = 'reflect'
+    elif mode == 'edge':
+        mode = 'nearest'
 
     def wrapped_func(arr):
         return function(arr, *extra_arguments, **extra_keywords)