diff --git a/lib/iris/analysis/__init__.py b/lib/iris/analysis/__init__.py
index 8ea451335d..122bc668f2 100644
--- a/lib/iris/analysis/__init__.py
+++ b/lib/iris/analysis/__init__.py
@@ -69,7 +69,7 @@
            'PEAK', 'PERCENTILE', 'PROPORTION', 'RMS', 'STD_DEV', 'SUM',
            'VARIANCE', 'WPERCENTILE', 'coord_comparison', 'Aggregator',
            'WeightedAggregator', 'clear_phenomenon_identity', 'Linear',
-           'AreaWeighted', 'Nearest')
+           'AreaWeighted', 'Nearest', 'UnstructuredNearest')
 
 
 class _CoordGroup(object):
@@ -2348,13 +2348,31 @@ def regridder(self, src_grid, target_grid):
 
 class UnstructuredNearest(object):
     """
-    This is a nearest-neighbour interpolation and regridding scheme for
-    regridding cubes whose latitude and longitude coordinates are mapped to the
-    same dimensions, rather than being orthogonal on independent dimensions.
+    This is a nearest-neighbour regridding scheme for regridding data whose
+    horizontal (X- and Y-axis) coordinates are mapped to the *same* dimensions,
+    rather than being orthogonal on independent dimensions.
 
-    Currently only supports regridding, not interpolation.
+    For latitude-longitude coordinates, the nearest-neighbour distances are
+    computed on the sphere, otherwise flat Euclidean distances are used.
+
+    The source X and Y coordinates can have any shape.
+
+    The target grid must be of the "normal" kind, i.e. it has separate,
+    1-dimensional X and Y coordinates.
+
+    Source and target XY coordinates must have the same coordinate system,
+    which may also be None.
+    If any of the XY coordinates are latitudes or longitudes, then they *all*
+    must be.  Otherwise, the corresponding X and Y coordinates must have the
+    same units in the source and grid cubes.
+
+    .. Note::
+        Currently only supports regridding, not interpolation.
 
     """
+    # Note: the argument requirements are simply those of the underlying
+    # regridder class,
+    # :class:`iris.analysis.trajectory.UnstructuredNearestNeigbourRegridder`.
     def __init__(self):
         """
         Nearest-neighbour interpolation and regridding scheme suitable for
@@ -2367,53 +2385,14 @@ def __init__(self):
     def __repr__(self):
         return 'UnstructuredNearest()'
 
-#    def interpolator(self, cube, coords):
-#        """
-#        Creates a nearest-neighbour interpolator to perform
-#        interpolation over the given :class:`~iris.cube.Cube` specified
-#        by the dimensions of the specified coordinates.
-#
-#        Typically you should use :meth:`iris.cube.Cube.interpolate` for
-#        interpolating a cube. There are, however, some situations when
-#        constructing your own interpolator is preferable. These are detailed
-#        in the :ref:`user guide <caching_an_interpolator>`.
-#
-#        Args:
-#
-#        * cube:
-#            The source :class:`iris.cube.Cube` to be interpolated.
-#        * coords:
-#            The names or coordinate instances that are to be
-#            interpolated over.
-#
-#        Returns:
-#            A callable with the interface:
-#
-#                `callable(sample_points, collapse_scalar=True)`
-#
-#            where `sample_points` is a sequence containing an array of values
-#            for each of the coordinates passed to this method, and
-#            `collapse_scalar` determines whether to remove length one
-#            dimensions in the result cube caused by scalar values in
-#            `sample_points`.
-#
-#            The values for coordinates that correspond to date/times
-#            may optionally be supplied as datetime.datetime or
-#            netcdftime.datetime instances.
-#
-#            For example, for the callable returned by:
-#            `Nearest().interpolator(cube, ['latitude', 'longitude'])`,
-#            sample_points must have the form
-#            `[new_lat_values, new_lon_values]`.
-#
-#        """
-#        return RectilinearInterpolator(cube, coords, 'nearest',
-#                                       self.extrapolation_mode)
+    # TODO: add interpolator usage
+    # def interpolator(self, cube):
 
     def regridder(self, src_cube, target_grid):
         """
-        Creates a nearest-neighbour regridder to perform regridding from the
-        source grid to the target grid.
+        Creates a nearest-neighbour regridder, of the
+        :class:`~iris.analysis.trajectory.UnstructuredNearestNeigbourRegridder`
+        type, to perform regridding from the source grid to the target grid.
 
         This can then be applied to any source data with the same structure as
         the original 'src_cube'.
@@ -2427,13 +2406,14 @@ def regridder(self, src_cube, target_grid):
 
         * src_cube:
             The :class:`~iris.cube.Cube` defining the source grid.
-            The X and Y coordinates must be mapped over the same dimensions.
+            The X and Y coordinates can have any shape, but must be mapped over
+            the same cube dimensions.
 
         * target_grid:
             The :class:`~iris.cube.Cube` defining the target grid.
-            It must have only 2 dimensions.
-            The X and Y coordinates must be one-dimensional and mapped to
-            different dimensions.
+            The X and Y coordinates must be one-dimensional dimension
+            coordinates, mapped to different dimensions.
+            All other cube components are ignored.
 
         Returns:
             A callable with the interface:
diff --git a/lib/iris/analysis/_interpolate_private.py b/lib/iris/analysis/_interpolate_private.py
index 69b11b8026..d25899c315 100644
--- a/lib/iris/analysis/_interpolate_private.py
+++ b/lib/iris/analysis/_interpolate_private.py
@@ -206,6 +206,11 @@ def _nearest_neighbour_indices_ndcoords(cube, sample_points, cache=None):
     Because this function can be slow for multidimensional coordinates,
     a 'cache' dictionary can be provided by the calling code.
 
+    .. Note::
+
+        If the points are longitudes/latitudes, these are handled correctly as
+        points on the sphere, but the values must be in 'degrees'.
+
     """
 
     # Developer notes:
diff --git a/lib/iris/analysis/trajectory.py b/lib/iris/analysis/trajectory.py
index 9ef7d7639f..fb5218100f 100644
--- a/lib/iris/analysis/trajectory.py
+++ b/lib/iris/analysis/trajectory.py
@@ -28,12 +28,14 @@
 
 import numpy as np
 
+from cf_units import Unit
 import iris.analysis
 import iris.coord_systems
 import iris.coords
 
 from iris.analysis._interpolate_private import \
     _nearest_neighbour_indices_ndcoords, linear as linear_regrid
+from iris.analysis._interpolation import snapshot_grid
 
 
 class _Segment(object):
@@ -278,24 +280,26 @@ def interpolate(cube, sample_points, method=None):
         fancy_source_indices = []
         region_slices = []
         n_index_length = len(column_indexes[0])
+        dims_reduced = [False] * n_index_length
         for i_ind in range(n_index_length):
             contents = [column_index[i_ind]
                         for column_index in column_indexes]
             each_used = [content != slice(None) for content in contents]
-            if not np.any(each_used):
-                # This dimension is not addressed by the operation.
-                # Use a ":" as the index.
-                fancy_index = slice(None)
-                # No sub-region selection for this dimension.
-                region_slice = slice(None)
-            elif np.all(each_used):
+            if np.all(each_used):
                 # This dimension is addressed : use a list of indices.
+                dims_reduced[i_ind] = True
                 # Select the region by min+max indices.
                 start_ind = np.min(contents)
                 stop_ind = 1 + np.max(contents)
                 region_slice = slice(start_ind, stop_ind)
                 # Record point indices with start subtracted from all of them.
                 fancy_index = list(np.array(contents) - start_ind)
+            elif not np.any(each_used):
+                # This dimension is not addressed by the operation.
+                # Use a ":" as the index.
+                fancy_index = slice(None)
+                # No sub-region selection for this dimension.
+                region_slice = slice(None)
             else:
                 # Should really never happen, if _ndcoords is right.
                 msg = ('Internal error in trajectory interpolation : point '
@@ -305,17 +309,34 @@ def interpolate(cube, sample_points, method=None):
             fancy_source_indices.append(fancy_index)
             region_slices.append(region_slice)
 
-        # NOTE: fetch the required (square-section) region of the source data.
-        # This is not quite as good as only fetching the individual points
-        # which are used, but it avoids creating a sub-cube for each point,
+        # Fetch the required (square-section) region of the source data.
+        # NOTE: This is not quite as good as only fetching the individual
+        # points used, but it avoids creating a sub-cube for each point,
         # which is very slow, especially when points are re-used a lot ...
         source_area_indices = tuple(region_slices)
         source_data = cube[source_area_indices].data
 
-        # Apply fancy indexing to get all the result data points.
+        # Transpose source data before indexing it to get the final result.
+        # Because.. the fancy indexing will replace the indexed (horizontal)
+        # dimensions with a new single dimension over trajectory points.
+        # Move those dimensions to the end *first* : this ensures that the new
+        # dimension also appears at the end, which is where we want it.
+        # Make a list of dims with the reduced ones last.
+        dims_reduced = np.array(dims_reduced)
+        dims_order = np.arange(n_index_length)
+        dims_order = np.concatenate((dims_order[~dims_reduced],
+                                     dims_order[dims_reduced]))
+        # Rearrange the data dimensions and the fancy indices into that order.
+        source_data = source_data.transpose(dims_order)
+        fancy_source_indices = [fancy_source_indices[i_dim]
+                                for i_dim in dims_order]
+
+        # Apply the fancy indexing to get all the result data points.
         source_data = source_data[fancy_source_indices]
+
         # "Fix" problems with missing datapoints producing odd values
         # when copied from a masked into an unmasked array.
+        # TODO: proper masked data handling.
         if np.ma.isMaskedArray(source_data):
             # This is **not** proper mask handling, because we cannot produce a
             # masked result, but it ensures we use a "filled" version of the
@@ -364,11 +385,13 @@ class UnstructuredNearestNeigbourRegridder(object):
     Encapsulate the operation of :meth:`iris.analysis.trajectory.interpolate`
     with given source and target grids.
 
-    TODO: cache the necessary bits of the operation so re-use can actually
-    be more efficient.
+    This is the type used by the :class:`~iris.analysis.UnstructuredNearest`
+    regridding scheme.
 
     """
-    def __init__(self, src_cube, target_grid):
+    # TODO: cache the necessary bits of the operation so re-use can actually
+    # be more efficient.
+    def __init__(self, src_cube, target_grid_cube):
         """
         A nearest-neighbour regridder to perform regridding from the source
         grid to the target grid.
@@ -380,13 +403,15 @@ def __init__(self, src_cube, target_grid):
 
         * src_cube:
             The :class:`~iris.cube.Cube` defining the source grid.
-            The X and Y coordinates must be mapped over the same dimensions.
+            The X and Y coordinates can have any shape, but must be mapped over
+            the same cube dimensions.
 
-        * target_grid:
-            The :class:`~iris.cube.Cube` defining the target grid.
-            It must have only 2 dimensions.
-            The X and Y coordinates must be one-dimensional and mapped to
-            different dimensions.
+        * target_grid_cube:
+            A :class:`~iris.cube.Cube`, whose X and Y coordinates specify a
+            desired target grid.
+            The X and Y coordinates must be one-dimensional dimension
+            coordinates, mapped to different dimensions.
+            All other cube components are ignored.
 
         Returns:
             regridder : (object)
@@ -395,78 +420,142 @@ def __init__(self, src_cube, target_grid):
                 `result_cube = regridder(data)`
 
             where `data` is a cube with the same grid as the original
-            `src_cube`, that is to be regridded to the `target_grid`.
+            `src_cube`, that is to be regridded to the `target_grid_cube`.
+
+        .. Note::
+
+            For latitude-longitude coordinates, the nearest-neighbour distances
+            are computed on the sphere, otherwise flat Euclidean distances are
+            used.
+
+            The source and target X and Y coordinates must all have the same
+            coordinate system, which may also be None.
+            If any X and Y coordinates are latitudes or longitudes, they *all*
+            must be.  Otherwise, the corresponding X and Y coordinates must
+            have the same units in the source and grid cubes.
 
         """
-        # Store the essential stuff
-        self.src_cube = src_cube
-        self.grid_cube = target_grid
-
-        # Quickly check the source data structure.
-        # TODO: replace asserts with code to raise user-intelligible errors.
-
-        # Has unique X and Y coords.
-        x_co = src_cube.coord(axis='x')
-        y_co = src_cube.coord(axis='y')
-        # They have a single common dimension, WHICH IS THE LAST.
-        src_ndim = src_cube.ndim
-        assert src_cube.coord_dims(x_co) == (src_ndim - 1,)
-        assert src_cube.coord_dims(y_co) == (src_ndim - 1,)
-
-        # Quickly check the target grid structure.
-        # TODO: ensure any errors are intelligible to the user.
-        # Has only 2 dims.
-        assert target_grid.ndim == 2
-        # Has unique X and Y coords.
-        x_co = target_grid.coord(axis='x')
-        y_co = target_grid.coord(axis='y')
-        # Each has a dimension to itself.
-        x_dims = target_grid.coord_dims(x_co)
-        y_dims = target_grid.coord_dims(y_co)
-        assert len(x_dims) == 1
-        assert len(y_dims) == 1
-        assert x_dims != y_dims
-
-        # Pre-calculate the sample points that will be needed.
-        # These are cast as a 'trajectory' to suit the method used.
-        x_vals = target_grid.coord('longitude').points
-        y_vals = target_grid.coord('latitude').points
-        x_2d, y_2d = np.meshgrid(x_vals, y_vals)
-        self.trajectory = (('longitude', x_2d.flatten()),
-                           ('latitude', y_2d.flatten()))
+        # Make a copy of the source cube, so we can convert coordinate units.
+        src_cube = src_cube.copy()
+
+        # Snapshot the target grid and check it is a "normal" grid.
+        tgt_x_coord, tgt_y_coord = snapshot_grid(target_grid_cube)
+
+        # Check that the source has unique X and Y coords over common dims.
+        if (not src_cube.coords(axis='x') or not src_cube.coords(axis='y')):
+            msg = 'Source cube must have X- and Y-axis coordinates.'
+            raise ValueError(msg)
+        src_x_coord = src_cube.coord(axis='x')
+        src_y_coord = src_cube.coord(axis='y')
+        if (src_cube.coord_dims(src_x_coord) !=
+                src_cube.coord_dims(src_y_coord)):
+            msg = ('Source cube X and Y coordinates must have the same '
+                   'cube dimensions.')
+            raise ValueError(msg)
+
+        # Record *copies* of the original grid coords, in the desired
+        # dimension order.
+        # This lets us convert the actual ones in use to units of "degrees".
+        self.src_grid_coords = [src_y_coord.copy(), src_x_coord.copy()]
+        self.tgt_grid_coords = [tgt_y_coord.copy(), tgt_x_coord.copy()]
+
+        # Check that all XY coords have suitable coordinate systems and units.
+        coords_all = [src_x_coord, src_y_coord, tgt_x_coord, tgt_y_coord]
+        cs = coords_all[0].coord_system
+        if not all(coord.coord_system == cs for coord in coords_all):
+            msg = ('Source and target cube X and Y coordinates must all have '
+                   'the same coordinate system.')
+            raise ValueError(msg)
+
+        # Check *all* X and Y coords are lats+lons, if any are.
+        latlons = ['latitude' in coord.name() or 'longitude' in coord.name()
+                   for coord in coords_all]
+        if any(latlons) and not all(latlons):
+            msg = ('If any X and Y coordinates are latitudes/longitudes, '
+                   'then they all must be.')
+            raise ValueError(msg)
+
+        self.grid_is_latlon = any(latlons)
+        if self.grid_is_latlon:
+            # Convert all XY coordinates to units of "degrees".
+            # N.B. already copied the target grid, so the result matches that.
+            for coord in coords_all:
+                try:
+                    coord.convert_units('degrees')
+                except ValueError:
+                    msg = ('Coordinate {!r} has units of {!r}, which does not '
+                           'convert to "degrees".')
+                    raise ValueError(msg.format(coord.name(),
+                                                str(coord.units)))
+        else:
+            # Check that source and target have the same X and Y units.
+            if (src_x_coord.units != tgt_x_coord.units or
+                    src_y_coord.units != tgt_y_coord.units):
+                msg = ('Source and target cube X and Y coordinates must '
+                       'have the same units.')
+                raise ValueError(msg)
+
+        # Record the resulting grid shape.
+        self.tgt_grid_shape = tgt_y_coord.shape + tgt_x_coord.shape
+
+        # Calculate sample points as 2d arrays, like broadcast (NY,1)*(1,NX).
+        x_2d, y_2d = np.meshgrid(tgt_x_coord.points, tgt_y_coord.points)
+        # Cast as a "trajectory", to suit the method used.
+        self.trajectory = ((tgt_x_coord.name(), x_2d.flatten()),
+                           (tgt_y_coord.name(), y_2d.flatten()))
 
     def __call__(self, src_cube):
-        # Check source cube matches original.
-        # For now, just a shape match will do.
-        # TODO: implement a more intelligent equivalence check.
-        # TODO: replace asserts with code to raise user-intelligible errors.
-        assert src_cube.shape == self.src_cube.shape
+        # Check the source cube X and Y coords match the original.
+        # Note: for now, this is sufficient to ensure a valid trajectory
+        # interpolation, but if in future we save + re-use the cache context
+        # for the 'interpolate' call, we may need more checks here.
+
+        # Check the given cube against the original.
+        x_cos = src_cube.coords(axis='x')
+        y_cos = src_cube.coords(axis='y')
+        if (not x_cos or not y_cos or
+                y_cos != [self.src_grid_coords[0]] or
+                x_cos != [self.src_grid_coords[1]]):
+            msg = ('The given cube is not defined on the same source '
+                   'grid as this regridder.')
+            raise ValueError(msg)
+
+        # Convert source XY coordinates to degrees if required.
+        if self.grid_is_latlon:
+            src_cube = src_cube.copy()
+            src_cube.coord(axis='x').convert_units('degrees')
+            src_cube.coord(axis='y').convert_units('degrees')
 
         # Get the basic interpolated results.
         result_trajectory_cube = interpolate(src_cube, self.trajectory,
                                              method='nearest')
 
         # Reconstruct this as a cube "like" the source data.
-        # TODO: sort out aux-coords, cell methods, cell measures ??
+        # TODO: handle all aux-coords, cell measures ??
 
-        # The shape is that of source data, minus the last dim, plus the target
-        # grid dimensions.
-        target_shape = (list(src_cube.shape)[:-1] + list(self.grid_cube.shape))
+        # The shape is that of the basic result, minus the trajectory (last)
+        # dimension, plus the target grid dimensions.
+        target_shape = result_trajectory_cube.shape[:-1] + self.tgt_grid_shape
         data_2d_x_and_y = result_trajectory_cube.data.reshape(target_shape)
 
         # Make a new result cube with the reshaped data.
         result_cube = iris.cube.Cube(data_2d_x_and_y)
         result_cube.metadata = src_cube.metadata
 
-        # Copy the 'preceding' dim coords from the source cube.
-        n_other_dims = src_cube.ndim - 1
-        for i_dim in range(n_other_dims):
-            co = src_cube.coord(dimensions=(i_dim,), dim_coords=True)
-            result_cube.add_dim_coord(co.copy(), i_dim)
-
-        # Copy the 'trailing' lat+lon coords from the grid cube.
-        for i_dim in (0, 1):
-            co = self.grid_cube.coord(dimensions=(i_dim,))
-            result_cube.add_dim_coord(co.copy(), i_dim + n_other_dims)
+        # Copy all the coords from the trajectory result.
+        i_trajectory_dim = result_trajectory_cube.ndim - 1
+        for coord in result_trajectory_cube.dim_coords:
+            dims = result_trajectory_cube.coord_dims(coord)
+            if i_trajectory_dim not in dims:
+                result_cube.add_dim_coord(coord.copy(), dims)
+        for coord in result_trajectory_cube.aux_coords:
+            dims = result_trajectory_cube.coord_dims(coord)
+            if i_trajectory_dim not in dims:
+                result_cube.add_aux_coord(coord.copy(), dims)
+
+        # Add the X+Y grid coords from the grid cube, mapped to the new Y and X
+        # dimensions, i.e. the last 2.
+        for i_dim, coord in enumerate(self.tgt_grid_coords):
+            result_cube.add_dim_coord(coord.copy(), i_dim + i_trajectory_dim)
 
         return result_cube
diff --git a/lib/iris/tests/integration/test_regridding.py b/lib/iris/tests/integration/test_regridding.py
index f12e79a8be..8255355ded 100644
--- a/lib/iris/tests/integration/test_regridding.py
+++ b/lib/iris/tests/integration/test_regridding.py
@@ -106,7 +106,7 @@ def setUp(self):
 
     def test_nearest(self):
         res = self.src.regrid(self.grid, UnstructuredNearest())
-        self.assertArrayShapeStats(res, (1, 6, 3, 4), 315.888737, 11.000729)
+        self.assertArrayShapeStats(res, (1, 6, 3, 4), 315.890808, 11.000724)
 
 
 if __name__ == "__main__":
diff --git a/lib/iris/tests/unit/analysis/trajectory/test_UnstructuredNearestNeighbourRegridder.py b/lib/iris/tests/unit/analysis/trajectory/test_UnstructuredNearestNeighbourRegridder.py
new file mode 100644
index 0000000000..c911e70bc2
--- /dev/null
+++ b/lib/iris/tests/unit/analysis/trajectory/test_UnstructuredNearestNeighbourRegridder.py
@@ -0,0 +1,306 @@
+# (C) British Crown Copyright 2016, Met Office
+#
+# This file is part of Iris.
+#
+# Iris is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by the
+# Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# Iris is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with Iris.  If not, see <http://www.gnu.org/licenses/>.
+"""
+Unit tests for
+:class:`iris.analysis.trajectory.UnstructuredNearestNeigbourRegridder`.
+
+"""
+from __future__ import (absolute_import, division, print_function)
+from six.moves import (filter, input, map, range, zip)  # noqa
+
+# Import iris.tests first so that some things can be initialised before
+# importing anything else.
+import iris.tests as tests
+
+from contextlib import contextmanager
+
+import numpy as np
+
+from iris.coords import AuxCoord, DimCoord
+from iris.coord_systems import GeogCS, RotatedGeogCS
+from iris.cube import Cube, CubeList
+import iris.tests.stock
+
+from iris.analysis.trajectory import \
+    UnstructuredNearestNeigbourRegridder as unn_gridder
+
+
+class MixinExampleSetup(object):
+    # Common code for regridder test classes.
+
+    def setUp(self):
+        # Basic test values.
+        src_x_y_value = np.array([
+            [20.12, 11.73, 0.01],
+            [120.23, -20.73, 1.12],
+            [290.34, 33.88, 2.23],
+            [-310.45, 57.8, 3.34]])
+        tgt_grid_x = np.array([-173.2, -100.3, -32.5, 1.4, 46.6, 150.7])
+        tgt_grid_y = np.array([-80.1, -30.2, 0.3, 47.4, 75.5])
+
+        # Make sample 1-D source cube.
+        src = Cube(src_x_y_value[:, 2])
+        src.add_aux_coord(AuxCoord(src_x_y_value[:, 0],
+                                   standard_name='longitude', units='degrees'),
+                          0)
+        src.add_aux_coord(AuxCoord(src_x_y_value[:, 1],
+                                   standard_name='latitude', units='degrees'),
+                          0)
+        self.src_cube = src
+
+        # Make sample grid cube.
+        grid = Cube(np.zeros(tgt_grid_y.shape + tgt_grid_x.shape))
+        grid.add_dim_coord(DimCoord(tgt_grid_y,
+                                    standard_name='latitude', units='degrees'),
+                           0)
+        grid.add_dim_coord(DimCoord(tgt_grid_x,
+                                    standard_name='longitude',
+                                    units='degrees'),
+                           1)
+        self.grid_cube = grid
+
+        # Make expected-result, from the expected source-index at each point.
+        expected_result_indices = np.array([
+            [1, 1, 1, 1, 1, 1],
+            [1, 2, 0, 0, 0, 1],
+            [1, 2, 2, 0, 0, 1],
+            [3, 2, 2, 3, 3, 3],
+            [3, 2, 3, 3, 3, 3]])
+        self.expected_data = self.src_cube.data[expected_result_indices]
+
+        # Make a 3D source cube, based on the existing 2d test data.
+        z_cubes = [src.copy() for _ in range(3)]
+        for i_z, z_cube in enumerate(z_cubes):
+            z_cube.add_aux_coord(DimCoord([i_z], long_name='z'))
+            z_cube.data = z_cube.data + 100.0 * i_z
+        self.src_z_cube = CubeList(z_cubes).merge_cube()
+
+        # Make a corresponding 3d expected result.
+        self.expected_data_zxy = \
+            self.src_z_cube.data[:, expected_result_indices]
+
+    def _check_expected(self, src_cube=None, grid_cube=None,
+                        expected_data=None,
+                        expected_coord_names=None):
+        # Test regridder creation + operation against expected results.
+        if src_cube is None:
+            src_cube = self.src_cube
+        if grid_cube is None:
+            grid_cube = self.grid_cube
+        gridder = unn_gridder(src_cube, grid_cube)
+        result = gridder(src_cube)
+        if expected_coord_names is not None:
+            # Check result coordinate identities.
+            self.assertEqual([coord.name() for coord in result.coords()],
+                             expected_coord_names)
+        if expected_data is None:
+            # By default, check against the 'standard' data result.
+            expected_data = self.expected_data
+        self.assertArrayEqual(result.data, expected_data)
+        return result
+
+
+class Test__init__(MixinExampleSetup, tests.IrisTest):
+    # Exercise all the constructor argument checks.
+
+    def test_fail_no_src_x(self):
+        self.src_cube.remove_coord('longitude')
+        msg_re = 'Source cube must have X- and Y-axis coordinates'
+        with self.assertRaisesRegexp(ValueError, msg_re):
+            unn_gridder(self.src_cube, self.grid_cube)
+
+    def test_fail_no_src_y(self):
+        self.src_cube.remove_coord('latitude')
+        msg_re = 'Source cube must have X- and Y-axis coordinates'
+        with self.assertRaisesRegexp(ValueError, msg_re):
+            unn_gridder(self.src_cube, self.grid_cube)
+
+    def test_fail_bad_src_dims(self):
+        self.src_cube = self.grid_cube
+        msg_re = 'Source.*same cube dimensions'
+        with self.assertRaisesRegexp(ValueError, msg_re):
+            unn_gridder(self.src_cube, self.grid_cube)
+
+    def test_fail_mixed_latlons(self):
+        self.src_cube.coord('longitude').rename('projection_x_coordinate')
+        msg_re = 'any.*latitudes/longitudes.*all must be'
+        with self.assertRaisesRegexp(ValueError, msg_re):
+            unn_gridder(self.src_cube, self.grid_cube)
+
+    def test_fail_bad_latlon_units(self):
+        self.grid_cube.coord('longitude').units = 'm'
+        msg_re = 'does not convert to "degrees"'
+        with self.assertRaisesRegexp(ValueError, msg_re):
+            unn_gridder(self.src_cube, self.grid_cube)
+
+    def test_fail_non_latlon_units_mismatch(self):
+        # Convert all to non-latlon system (does work: see in "Test__call__").
+        for cube in (self.src_cube, self.grid_cube):
+            for axis_name in ('x', 'y'):
+                coord = cube.coord(axis=axis_name)
+                coord_name = 'projection_{}_coordinate'.format(axis_name)
+                coord.rename(coord_name)
+                coord.units = 'm'
+        # Change one of the output units.
+        self.grid_cube.coord(axis='x').units = '1'
+        msg_re = 'Source and target.*must have the same units'
+        with self.assertRaisesRegexp(ValueError, msg_re):
+            unn_gridder(self.src_cube, self.grid_cube)
+
+    def test_fail_no_tgt_x(self):
+        self.grid_cube.remove_coord('longitude')
+        msg_re = 'must contain a single 1D x coordinate'
+        with self.assertRaisesRegexp(ValueError, msg_re):
+            unn_gridder(self.src_cube, self.grid_cube)
+
+    def test_fail_no_tgt_y(self):
+        self.grid_cube.remove_coord('latitude')
+        msg_re = 'must contain a single 1D y coordinate'
+        with self.assertRaisesRegexp(ValueError, msg_re):
+            unn_gridder(self.src_cube, self.grid_cube)
+
+    def test_fail_src_cs_mismatch(self):
+        cs = GeogCS(1000.0)
+        self.src_cube.coord('latitude').coord_system = cs
+        msg_re = 'must all have the same coordinate system'
+        with self.assertRaisesRegexp(ValueError, msg_re):
+            unn_gridder(self.src_cube, self.grid_cube)
+
+    def test_fail_tgt_cs_mismatch(self):
+        cs = GeogCS(1000.0)
+        self.grid_cube.coord('latitude').coord_system = cs
+        msg_re = 'x.*and y.*must have the same coordinate system'
+        with self.assertRaisesRegexp(ValueError, msg_re):
+            unn_gridder(self.src_cube, self.grid_cube)
+
+    def test_fail_src_tgt_cs_mismatch(self):
+        cs = GeogCS(1000.0)
+        self.src_cube.coord('latitude').coord_system = cs
+        self.src_cube.coord('longitude').coord_system = cs
+        msg_re = 'Source and target.*same coordinate system'
+        with self.assertRaisesRegexp(ValueError, msg_re):
+            unn_gridder(self.src_cube, self.grid_cube)
+
+
+class Test__call__(MixinExampleSetup, tests.IrisTest):
+    # Test regridder operation and results.
+
+    def test_basic_latlon(self):
+        # Check a test operation.
+        self._check_expected(expected_coord_names=['latitude', 'longitude'],
+                             expected_data=self.expected_data)
+
+    def test_non_latlon(self):
+        # Check different answer in cartesian coordinates (no wrapping, etc).
+        # Convert to non-latlon system, with the same coord values.
+        for cube in (self.src_cube, self.grid_cube):
+            for axis_name in ('x', 'y'):
+                coord = cube.coord(axis=axis_name)
+                coord_name = 'projection_{}_coordinate'.format(axis_name)
+                coord.rename(coord_name)
+                coord.units = 'm'
+        # Check for a somewhat different result.
+        non_latlon_indices = np.array([
+            [3, 0, 0, 0, 1, 1],
+            [3, 0, 0, 0, 0, 1],
+            [3, 0, 0, 0, 0, 1],
+            [3, 0, 0, 0, 0, 1],
+            [3, 0, 0, 0, 0, 1]])
+        expected_data = self.src_cube.data[non_latlon_indices]
+        self._check_expected(expected_data=expected_data)
+
+    def test_multidimensional_xy(self):
+        # Recast the 4-point source cube as 2*2 : should yield the same result.
+        co_x = self.src_cube.coord(axis='x')
+        co_y = self.src_cube.coord(axis='y')
+        new_src = Cube(self.src_cube.data.reshape((2, 2)))
+        new_x_co = AuxCoord(co_x.points.reshape((2, 2)),
+                            standard_name='longitude', units='degrees')
+        new_y_co = AuxCoord(co_y.points.reshape((2, 2)),
+                            standard_name='latitude', units='degrees')
+        new_src.add_aux_coord(new_x_co, (0, 1))
+        new_src.add_aux_coord(new_y_co, (0, 1))
+        self._check_expected(src_cube=new_src)
+
+    def test_transposed_grid(self):
+        # Show that changing the order of the grid X and Y has no effect.
+        new_grid_cube = self.grid_cube.copy()
+        new_grid_cube.transpose((1, 0))
+        # Check that the new grid is in (X, Y) order.
+        self.assertEqual([coord.name() for coord in new_grid_cube.coords()],
+                         ['longitude', 'latitude'])
+        # Check that the result is the same, dimension order is still Y,X.
+        self._check_expected(grid_cube=new_grid_cube,
+                             expected_coord_names=['latitude', 'longitude'])
+
+    def test_compatible_source(self):
+        # Check operation on data with different dimensions to the original
+        # source cube for the regridder creation.
+        gridder = unn_gridder(self.src_cube, self.grid_cube)
+        result = gridder(self.src_z_cube)
+        self.assertEqual([coord.name() for coord in result.coords()],
+                         ['z', 'latitude', 'longitude'])
+        self.assertArrayEqual(result.data, self.expected_data_zxy)
+
+    def test_fail_incompatible_source(self):
+        # Check that a slightly modified source cube is *not* acceptable.
+        modified_src_cube = self.src_cube.copy()
+        points = modified_src_cube.coord(axis='x').points
+        points[0] += 0.01
+        modified_src_cube.coord(axis='x').points = points
+        gridder = unn_gridder(self.src_cube, self.grid_cube)
+        msg = 'not defined on the same source grid'
+        with self.assertRaisesRegexp(ValueError, msg):
+            gridder(modified_src_cube)
+
+    def test_transposed_source(self):
+        # Check operation on data where the 'trajectory' dimension is not the
+        # last one.
+        src_z_cube = self.src_z_cube
+        src_z_cube.transpose((1, 0))
+        self._check_expected(src_cube=src_z_cube,
+                             expected_data=self.expected_data_zxy)
+
+    def test_radians_degrees(self):
+        # Check source + target unit conversions, grid and result in degrees.
+        for axis_name in ('x', 'y'):
+            self.src_cube.coord(axis=axis_name).convert_units('radians')
+            self.grid_cube.coord(axis=axis_name).convert_units('degrees')
+        result = self._check_expected()
+        self.assertEqual(result.coord(axis='x').units, 'degrees')
+
+    def test_degrees_radians(self):
+        # Check source + target unit conversions, grid and result in radians.
+        for axis_name in ('x', 'y'):
+            self.src_cube.coord(axis=axis_name).convert_units('degrees')
+            self.grid_cube.coord(axis=axis_name).convert_units('radians')
+        result = self._check_expected()
+        self.assertEqual(result.coord(axis='x').units, 'radians')
+
+    def test_alternative_cs(self):
+        # Check the result is just the same in a different coordinate system.
+        cs = RotatedGeogCS(grid_north_pole_latitude=75.3,
+                           grid_north_pole_longitude=102.5,
+                           ellipsoid=GeogCS(100.0))
+        for cube in (self.src_cube, self.grid_cube):
+            for coord_name in ('longitude', 'latitude'):
+                cube.coord(coord_name).coord_system = cs
+        self._check_expected()
+
+
+if __name__ == "__main__":
+    tests.main()