diff --git a/lib/iris/analysis/cartography.py b/lib/iris/analysis/cartography.py
index 1f98e2377d..a7447907ba 100644
--- a/lib/iris/analysis/cartography.py
+++ b/lib/iris/analysis/cartography.py
@@ -713,3 +713,379 @@ def project(cube, target_proj, nx=None, ny=None):
     new_cube.metadata = cube.metadata
 
     return new_cube, extent
+
+
+def _transform_xy(crs_from, x, y, crs_to):
+    """
+    Shorthand function to transform 2d points between coordinate
+    reference systems.
+
+    Args:
+
+    * crs_from, crs_to (:class:`cartopy.crs.Projection`):
+        The coordinate reference systems.
+    * x, y (arrays):
+        point locations defined in 'crs_from'.
+
+    Returns:
+        x, y :  Arrays of locations defined in 'crs_to'.
+
+    """
+    pts = crs_to.transform_points(crs_from, x, y)
+    return pts[..., 0], pts[..., 1]
+
+
+def _inter_crs_differentials(crs1, x, y, crs2):
+    """
+    Calculate coordinate partial differentials from crs1 to crs2.
+
+    Returns dx2/dx1, dy2/dx1, dx2/dy1 and dy2/dy1, at given locations.
+
+    Args:
+
+    * crs1, crs2 (`cartopy.crs.Projection`):
+        The coordinate systems, "from" and "to".
+    * x, y (array):
+        Point locations defined in 'crs1'.
+
+    Returns:
+        (dx2/dx1, dy2/dx1, dx2/dy1, dy2/dy1) at given locations. Each
+        element of this tuple will be the same shape as the 'x' and 'y'
+        arrays and will be the partial differentials between the two systems.
+
+    """
+    # Get locations in target crs.
+    crs2_x, crs2_y = _transform_xy(crs1, x, y, crs2)
+
+    # Define small x-deltas in the source crs.
+    VECTOR_DELTAS_FACTOR = 360000.0  # Empirical factor to obtain small delta.
+    delta_x = (crs1.x_limits[1] - crs1.x_limits[0]) / VECTOR_DELTAS_FACTOR
+    delta_x = delta_x * np.ones(x.shape)
+    eps = 1e-9
+    # Reverse deltas where we would otherwise step outside the valid range.
+    invalid_dx = x + delta_x > crs1.x_limits[1] - eps
+    delta_x[invalid_dx] = -delta_x[invalid_dx]
+    # Calculate the transformed point with x = x + dx.
+    crs2_x2, crs2_y2 = _transform_xy(crs1, x + delta_x, y, crs2)
+    # Form differentials wrt dx.
+    dx2_dx = (crs2_x2 - crs2_x) / delta_x
+    dy2_dx = (crs2_y2 - crs2_y) / delta_x
+
+    # Define small y-deltas in the source crs.
+    delta_y = (crs1.y_limits[1] - crs1.y_limits[0]) / VECTOR_DELTAS_FACTOR
+    delta_y = delta_y * np.ones(y.shape)
+    # Reverse deltas where we would otherwise step outside the valid range.
+    invalid_dy = y + delta_y > crs1.y_limits[1] - eps
+    delta_y[invalid_dy] = -delta_y[invalid_dy]
+    # Calculate the transformed point with y = y + dy.
+    crs2_x2, crs2_y2 = _transform_xy(crs1, x, y + delta_y, crs2)
+    # Form differentials wrt dy.
+    dx2_dy = (crs2_x2 - crs2_x) / delta_y
+    dy2_dy = (crs2_y2 - crs2_y) / delta_y
+
+    return dx2_dx, dy2_dx, dx2_dy, dy2_dy
+
+
+def _crs_distance_differentials(crs, x, y):
+    """
+    Calculate d(distance) / d(x) and ... / d(y) for a coordinate
+    reference system at specified locations.
+
+    Args:
+
+    * crs (:class:`cartopy.crs.Projection`):
+        The coordinate reference system.
+    * x, y (array):
+        Locations at which to calculate the differentials,
+        defined in 'crs' coordinate reference system.
+
+    Returns:
+        (abs(ds/dx), abs(ds/dy)).
+        Numerically approximated partial differentials,
+        i.e. scaling factors between changes in distance and changes in
+        coordinate values.
+
+    """
+    # Make a true-latlon coordinate system for distance calculations.
+    crs_latlon = ccrs.Geodetic(globe=ccrs.Globe(ellipse='sphere'))
+    # Transform points to true-latlon (just to get the true latitudes).
+    _, true_lat = _transform_xy(crs, x, y, crs_latlon)
+    # Get coordinate differentials w.r.t. true-latlon.
+    dlon_dx, dlat_dx, dlon_dy, dlat_dy = \
+        _inter_crs_differentials(crs, x, y, crs_latlon)
+    # Calculate effective scalings of X and Y coordinates.
+    lat_factor = np.cos(np.deg2rad(true_lat))**2
+    ds_dx = np.sqrt(dlat_dx * dlat_dx + dlon_dx * dlon_dx * lat_factor)
+    ds_dy = np.sqrt(dlat_dy * dlat_dy + dlon_dy * dlon_dy * lat_factor)
+    return ds_dx, ds_dy
+
+
+def _transform_distance_vectors(src_crs, x, y, u_dist, v_dist, tgt_crs):
+    """
+    Transform distance vectors from one coordinate reference system to
+    another, preserving magnitude and physical direction.
+
+    Args:
+
+    * src_crs, tgt_crs (`cartopy.crs.Projection`):
+        The source and target coordinate reference systems.
+    * x, y (array):
+        Locations of each vector defined in 'src_crs'.
+    * u_dist, v_dist (array):
+        Components of each vector along the x and y directions of 'src_crs'
+        at each location.
+
+    Returns:
+        (ut_dist, vt_dist): Tuple of arrays containing the vector components
+        along the x and y directions of 'tgt_crs' at each location.
+
+    """
+    # Get distance scalings for source crs.
+    ds_dx1, ds_dy1 = _crs_distance_differentials(src_crs, x, y)
+    # Get distance scalings for target crs.
+    x2, y2 = _transform_xy(src_crs, x, y, tgt_crs)
+    ds_dx2, ds_dy2 = _crs_distance_differentials(tgt_crs, x2, y2)
+    # Scale input distance vectors --> source-coordinate differentials.
+    u1, v1 = u_dist / ds_dx1, v_dist / ds_dy1
+    # Transform vectors into the target system.
+    dx2_x1, dy2_x1, dx2_y1, dy2_y1 = _inter_crs_differentials(
+        src_crs, x, y, tgt_crs)
+    u2 = dx2_x1 * u1 + dx2_y1 * v1
+    v2 = dy2_x1 * u1 + dy2_y1 * v1
+    # Rescale output coordinate vectors --> target distance vectors.
+    u2_dist, v2_dist = u2 * ds_dx2, v2 * ds_dy2
+
+    return u2_dist, v2_dist
+
+
+def _transform_distance_vectors_tolerance_mask(src_crs, x, y, tgt_crs):
+    """
+    Return a mask that can be applied to data array to mask elements
+    where the magnitude of vectors are not preserved due to numerical
+    errors introduced by the tranformation between coordinate systems.
+
+    Args:
+    * src_crs (`cartopy.crs.Projection`):
+        The source coordinate reference systems.
+    * x, y (array):
+        Locations of each vector defined in 'src_crs'.
+    * tgt_crs (`cartopy.crs.Projection`):
+        The target coordinate reference systems.
+
+    Returns:
+        2d boolean array that is the same shape as x and y.
+
+    """
+    if x.shape != y.shape:
+        raise ValueError('Arrays do not have matching shapes. '
+                         'x.shape is {}, y.shape is {}.'.format(x.shape,
+                                                                y.shape))
+    ones = np.ones(x.shape)
+    zeros = np.zeros(x.shape)
+    u_one_t, v_zero_t = _transform_distance_vectors(src_crs, x, y,
+                                                    ones, zeros, tgt_crs)
+    u_zero_t, v_one_t = _transform_distance_vectors(src_crs, x, y,
+                                                    zeros, ones, tgt_crs)
+    # Squared magnitudes should be equal to one within acceptable tolerance.
+    # A value of atol=2e-3 is used, which corresponds to a change in magnitude
+    # of approximately 0.1%.
+    sqmag_1_0 = u_one_t**2 + v_zero_t**2
+    sqmag_0_1 = u_zero_t**2 + v_one_t**2
+    mask = np.logical_not(
+        np.logical_and(np.isclose(sqmag_1_0, ones, atol=2e-3),
+                       np.isclose(sqmag_0_1, ones, atol=2e-3)))
+    return mask
+
+
+def rotate_winds(u_cube, v_cube, target_cs):
+    """
+    Transform wind vectors to a different coordinate system.
+
+    The input cubes contain U and V components parallel to the local X and Y
+    directions of the input grid at each point.
+
+    The output cubes contain the same winds, at the same locations, but
+    relative to the grid directions of a different coordinate system.
+    Thus in vector terms, the magnitudes will always be the same, but the
+    angles can be different.
+
+    The outputs retain the original horizontal dimension coordinates, but
+    also have two 2-dimensional auxiliary coordinates containing the X and
+    Y locations in the target coordinate system.
+
+    Args:
+
+    * u_cube
+        An instance of :class:`iris.cube.Cube` that contains the x-component
+        of the vector.
+    * v_cube
+        An instance of :class:`iris.cube.Cube` that contains the y-component
+        of the vector.
+    * target_cs
+        An instance of :class:`iris.coord_systems.CoordSystem` that specifies
+        the new grid directions.
+
+    Returns:
+        A (u', v') tuple of :class:`iris.cube.Cube` instances that are the u
+        and v components in the requested target coordinate system.
+        The units are the same as the inputs.
+
+    .. note::
+
+        The U and V values relate to distance, with units such as 'm s-1'.
+        These are not the same as coordinate vectors, which transform in a
+        different manner.
+
+    .. note::
+
+        The names of the output cubes are those of the inputs, prefixed with
+        'transformed\_' (e.g. 'transformed_x_wind').
+
+    .. warning::
+
+        Conversion between rotated-pole and non-rotated systems can be
+        expressed analytically.  However, this function always uses a numerical
+        approach. In locations where this numerical approach does not preserve
+        magnitude to an accuracy of 0.1%, the corresponding elements of the
+        returned cubes will be masked.
+
+    """
+    # Check u_cube and v_cube have the same shape. We iterate through
+    # the u and v cube slices which relies on the shapes matching.
+    if u_cube.shape != v_cube.shape:
+        msg = 'Expected u and v cubes to have the same shape. ' \
+              'u cube has shape {}, v cube has shape {}.'
+        raise ValueError(msg.format(u_cube.shape, v_cube.shape))
+
+    # Check the u_cube and v_cube have the same x and y coords.
+    msg = 'Coordinates differ between u and v cubes. Coordinate {!r} from ' \
+          'u cube does not equal coordinate {!r} from v cube.'
+    if u_cube.coord(axis='x') != v_cube.coord(axis='x'):
+        raise ValueError(msg.format(u_cube.coord(axis='x').name(),
+                                    v_cube.coord(axis='x').name()))
+    if u_cube.coord(axis='y') != v_cube.coord(axis='y'):
+        raise ValueError(msg.format(u_cube.coord(axis='y').name(),
+                                    v_cube.coord(axis='y').name()))
+
+    # Check x and y coords have the same coordinate system.
+    x_coord = u_cube.coord(axis='x')
+    y_coord = u_cube.coord(axis='y')
+    if x_coord.coord_system != y_coord.coord_system:
+        msg = "Coordinate systems of x and y coordinates differ. " \
+              "Coordinate {!r} has a coord system of {!r}, but coordinate " \
+              "{!r} has a coord system of {!r}."
+        raise ValueError(msg.format(x_coord.name(), x_coord.coord_system,
+                                    y_coord.name(), y_coord.coord_system))
+
+    # Convert from iris coord systems to cartopy CRSs to access
+    # transform functionality. Use projection as cartopy
+    # transform_vectors relies on x_limits and y_limits.
+    if x_coord.coord_system is not None:
+        src_crs = x_coord.coord_system.as_cartopy_projection()
+    else:
+        # Default to Geodetic (but actually use PlateCarree as a
+        # projection is needed).
+        src_crs = ccrs.PlateCarree()
+    target_crs = target_cs.as_cartopy_projection()
+
+    # Check the number of dimensions of the x and y coords is the same.
+    # Subsequent logic assumes either both 1d or both 2d.
+    x = x_coord.points
+    y = y_coord.points
+    if x.ndim != y.ndim or x.ndim > 2 or y.ndim > 2:
+        msg = 'x and y coordinates must have the same number of dimensions ' \
+              'and be either 1D or 2D. The number of dimensions are {} and ' \
+              '{}, respectively.'.format(x.ndim, y.ndim)
+        raise ValueError(msg)
+
+    # Check the dimension mappings match between u_cube and v_cube.
+    if u_cube.coord_dims(x_coord) != v_cube.coord_dims(x_coord):
+        raise ValueError('Dimension mapping of x coordinate differs '
+                         'between u and v cubes.')
+    if u_cube.coord_dims(y_coord) != v_cube.coord_dims(y_coord):
+        raise ValueError('Dimension mapping of y coordinate differs '
+                         'between u and v cubes.')
+    x_dims = u_cube.coord_dims(x_coord)
+    y_dims = u_cube.coord_dims(y_coord)
+
+    # Convert points to 2D, if not already, and determine dims.
+    if x.ndim == y.ndim == 1:
+        x, y = np.meshgrid(x, y)
+        dims = (y_dims[0], x_dims[0])
+    else:
+        dims = x_dims
+
+    # Transpose x, y 2d arrays to match the order in cube's data
+    # array so that x, y and the sliced data all line up.
+    if dims[0] > dims[1]:
+        x = x.transpose()
+        y = y.transpose()
+
+    # Create resulting cubes.
+    ut_cube = u_cube.copy()
+    vt_cube = v_cube.copy()
+    ut_cube.rename('transformed_{}'.format(u_cube.name()))
+    vt_cube.rename('transformed_{}'.format(v_cube.name()))
+
+    # Calculate mask based on preservation of magnitude.
+    mask = _transform_distance_vectors_tolerance_mask(src_crs, x, y,
+                                                      target_crs)
+    apply_mask = mask.any()
+    if apply_mask:
+        # Make masked arrays to accept masking.
+        ut_cube.data = ma.asanyarray(ut_cube.data)
+        vt_cube.data = ma.asanyarray(vt_cube.data)
+
+    # Project vectors with u, v components one horiz slice at a time and
+    # insert into the resulting cubes.
+    shape = list(u_cube.shape)
+    for dim in dims:
+        shape[dim] = 1
+    ndindex = np.ndindex(*shape)
+    for index in ndindex:
+        index = list(index)
+        for dim in dims:
+            index[dim] = slice(None, None)
+        index = tuple(index)
+        u = u_cube.data[index]
+        v = v_cube.data[index]
+        ut, vt = _transform_distance_vectors(src_crs, x, y, u, v, target_crs)
+        if apply_mask:
+            ut = ma.asanyarray(ut)
+            ut[mask] = ma.masked
+            vt = ma.asanyarray(vt)
+            vt[mask] = ma.masked
+        ut_cube.data[index] = ut
+        vt_cube.data[index] = vt
+
+    # Calculate new coords of locations in target coordinate system.
+    xyz_tran = target_crs.transform_points(src_crs, x, y)
+    xt = xyz_tran[..., 0].reshape(x.shape)
+    yt = xyz_tran[..., 1].reshape(y.shape)
+
+    # Transpose xt, yt 2d arrays to match the dim order
+    # of the original x an y arrays - i.e. undo the earlier
+    # transpose (if applied).
+    if dims[0] > dims[1]:
+        xt = xt.transpose()
+        yt = yt.transpose()
+
+    xt_coord = iris.coords.AuxCoord(xt,
+                                    standard_name='projection_x_coordinate',
+                                    coord_system=target_cs)
+    yt_coord = iris.coords.AuxCoord(yt,
+                                    standard_name='projection_y_coordinate',
+                                    coord_system=target_cs)
+    # Set units based on coord_system.
+    if isinstance(target_cs, (iris.coord_systems.GeogCS,
+                              iris.coord_systems.RotatedGeogCS)):
+        xt_coord.units = yt_coord.units = 'degrees'
+    else:
+        xt_coord.units = yt_coord.units = 'm'
+
+    ut_cube.add_aux_coord(xt_coord, dims)
+    ut_cube.add_aux_coord(yt_coord, dims)
+    vt_cube.add_aux_coord(xt_coord.copy(), dims)
+    vt_cube.add_aux_coord(yt_coord.copy(), dims)
+
+    return ut_cube, vt_cube
diff --git a/lib/iris/coord_systems.py b/lib/iris/coord_systems.py
index af5bb63449..a66dbc930f 100644
--- a/lib/iris/coord_systems.py
+++ b/lib/iris/coord_systems.py
@@ -1,4 +1,4 @@
-# (C) British Crown Copyright 2010 - 2014, Met Office
+# (C) British Crown Copyright 2010 - 2015, Met Office
 #
 # This file is part of Iris.
 #
@@ -315,12 +315,20 @@ def xml_element(self, doc):
         return CoordSystem.xml_element(self, doc, self._pretty_attrs())
 
     def as_cartopy_crs(self):
+        globe = None
+        if self.ellipsoid is not None:
+            globe = self.ellipsoid.as_cartopy_globe()
         return ccrs.RotatedGeodetic(self.grid_north_pole_longitude,
-                                    self.grid_north_pole_latitude)
+                                    self.grid_north_pole_latitude,
+                                    globe=globe)
 
     def as_cartopy_projection(self):
+        globe = None
+        if self.ellipsoid is not None:
+            globe = self.ellipsoid.as_cartopy_globe()
         return ccrs.RotatedPole(self.grid_north_pole_longitude,
-                                self.grid_north_pole_latitude)
+                                self.grid_north_pole_latitude,
+                                globe=globe)
 
 
 class TransverseMercator(CoordSystem):
diff --git a/lib/iris/tests/test_coding_standards.py b/lib/iris/tests/test_coding_standards.py
index 16fdc6a399..dde4b7a6bf 100644
--- a/lib/iris/tests/test_coding_standards.py
+++ b/lib/iris/tests/test_coding_standards.py
@@ -1,4 +1,4 @@
-# (C) British Crown Copyright 2013 - 2014, Met Office
+# (C) British Crown Copyright 2013 - 2015, Met Office
 #
 # This file is part of Iris.
 #
@@ -114,7 +114,6 @@ class StandardReportWithExclusions(pep8.StandardReport):
         '*/iris/tests/test_io_init.py',
         '*/iris/tests/test_iterate.py',
         '*/iris/tests/test_load.py',
-        '*/iris/tests/test_mapping.py',
         '*/iris/tests/test_merge.py',
         '*/iris/tests/test_pickling.py',
         '*/iris/tests/test_pp_cf.py',
diff --git a/lib/iris/tests/test_mapping.py b/lib/iris/tests/test_mapping.py
index fde77e1e36..cb5092f129 100644
--- a/lib/iris/tests/test_mapping.py
+++ b/lib/iris/tests/test_mapping.py
@@ -18,10 +18,10 @@
 Tests map creation.
 
 """
-
 from __future__ import (absolute_import, division, print_function)
 
-# import iris tests first so that some things can be initialised before importing anything else
+# import iris tests first so that some things can be initialised before
+# importing anything else
 import iris.tests as tests
 
 import numpy as np
@@ -39,6 +39,12 @@
     import iris.plot as iplt
 
 
+# A specific cartopy Globe matching the iris RotatedGeogCS default.
+_DEFAULT_GLOBE = ccrs.Globe(semimajor_axis=6371229.0,
+                            semiminor_axis=6371229.0,
+                            ellipse=None)
+
+
 @tests.skip_plot
 class TestBasic(tests.GraphicsTest):
     cube = iris.tests.stock.realistic_4d()
@@ -61,13 +67,13 @@ def test_unmappable(self):
 
     def test_default_projection_and_extent(self):
         self.assertEqual(iplt.default_projection(self.cube),
-                         ccrs.RotatedPole(357.5 - 180, 37.5)
-                         )
+                         ccrs.RotatedPole(357.5 - 180, 37.5,
+                                          globe=_DEFAULT_GLOBE))
 
-        np_testing.assert_array_almost_equal(iplt.default_projection_extent(self.cube),
-                                             (3.59579163e+02, 3.59669159e+02, -1.28250003e-01, -3.82499993e-02),
-                                             decimal=3
-                                             )
+        np_testing.assert_array_almost_equal(
+            iplt.default_projection_extent(self.cube),
+            (3.59579163e+02, 3.59669159e+02, -1.28250003e-01, -3.82499993e-02),
+            decimal=3)
 
 
 @tests.skip_data
@@ -78,8 +84,14 @@ def setUp(self):
 
         # Make a cube that can't be located on the globe.
         cube = iris.cube.Cube(src_cube.data)
-        cube.add_dim_coord(iris.coords.DimCoord(np.arange(96, dtype=np.float32) * 100, long_name='x', units='m'), 1)
-        cube.add_dim_coord(iris.coords.DimCoord(np.arange(73, dtype=np.float32) * 100, long_name='y', units='m'), 0)
+        cube.add_dim_coord(
+            iris.coords.DimCoord(np.arange(96, dtype=np.float32) * 100,
+                                 long_name='x', units='m'),
+            1)
+        cube.add_dim_coord(
+            iris.coords.DimCoord(np.arange(73, dtype=np.float32) * 100,
+                                 long_name='y', units='m'),
+            0)
         cube.standard_name = 'air_temperature'
         cube.units = 'K'
         cube.assert_valid()
@@ -94,7 +106,8 @@ def test_simple(self):
 @tests.skip_plot
 class TestMappingSubRegion(tests.GraphicsTest):
     def setUp(self):
-        cube_path = tests.get_data_path(('PP', 'aPProt1', 'rotatedMHtimecube.pp'))
+        cube_path = tests.get_data_path(
+            ('PP', 'aPProt1', 'rotatedMHtimecube.pp'))
         cube = iris.load_cube(cube_path)[0]
         # make the data smaller to speed things up.
         self.cube = cube[::10, ::10]
@@ -129,12 +142,13 @@ def test_simple(self):
 
     def test_default_projection_and_extent(self):
         self.assertEqual(iplt.default_projection(self.cube),
-                          ccrs.RotatedPole(357.5 - 180, 37.5)
-                          )
+                         ccrs.RotatedPole(357.5 - 180, 37.5,
+                                          globe=_DEFAULT_GLOBE))
+
+        np_testing.assert_array_almost_equal(
+            iplt.default_projection_extent(self.cube),
+            (313.01998901, 391.11999512, -22.48999977, 24.80999947))
 
-        np_testing.assert_array_almost_equal(iplt.default_projection_extent(self.cube),
-                                             (313.01998901, 391.11999512, -22.48999977, 24.80999947)
-                                             )
 
 @tests.skip_data
 @tests.skip_plot
@@ -143,7 +157,8 @@ def setUp(self):
         self.cube = iris.tests.stock.global_pp()
         self.few = 4
         self.few_levels = range(280, 300, 5)
-        self.many_levels = np.linspace(self.cube.data.min(), self.cube.data.max(), 40)
+        self.many_levels = np.linspace(
+            self.cube.data.min(), self.cube.data.max(), 40)
 
     def test_simple(self):
         iplt.contour(self.cube)
@@ -172,13 +187,15 @@ def test_keywords(self):
 class TestBoundedCube(tests.GraphicsTest):
     def setUp(self):
         self.cube = iris.tests.stock.global_pp()
-        # Add some bounds to this data (this will actually make the bounds invalid as they
-        # will straddle the north pole and overlap on the date line, but that doesn't matter for this test.)
+        # Add some bounds to this data (this will actually make the bounds
+        # invalid as they will straddle the north pole and overlap on the
+        # dateline, but that doesn't matter for this test.)
         self.cube.coord('latitude').guess_bounds()
         self.cube.coord('longitude').guess_bounds()
 
     def test_pcolormesh(self):
-        # pcolormesh can only be drawn in native coordinates (or more specifically, in coordinates that don't wrap).
+        # pcolormesh can only be drawn in native coordinates (or more
+        # specifically, in coordinates that don't wrap).
         plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
         iplt.pcolormesh(self.cube)
         self.check_graphic()
@@ -191,12 +208,12 @@ def test_default_projection_and_extent(self):
         self.assertEqual(iplt.default_projection(self.cube),
                          ccrs.PlateCarree())
         np_testing.assert_array_almost_equal(
-             iplt.default_projection_extent(self.cube),
-             [0., 360., -89.99995422, 89.99998474])
+            iplt.default_projection_extent(self.cube),
+            [0., 360., -89.99995422, 89.99998474])
         np_testing.assert_array_almost_equal(
-             iplt.default_projection_extent(
-                 self.cube, mode=iris.coords.BOUND_MODE),
-                 [-1.875046, 358.124954, -91.24995422, 91.24998474])
+            iplt.default_projection_extent(
+                self.cube, mode=iris.coords.BOUND_MODE),
+            [-1.875046, 358.124954, -91.24995422, 91.24998474])
 
 
 @tests.skip_data
@@ -214,8 +231,8 @@ def test_pcolormesh(self):
 
     def test_grid(self):
         iplt.pcolormesh(self.cube, facecolors='none', edgecolors='blue')
-        # the result is a graphic which has coloured edges. This is a mpl bug, see
-        # https://github.com/matplotlib/matplotlib/issues/1302
+        # the result is a graphic which has coloured edges. This is a mpl bug,
+        # see https://github.com/matplotlib/matplotlib/issues/1302
         self.check_graphic()
 
     def test_outline(self):
diff --git a/lib/iris/tests/unit/analysis/cartography/test_rotate_winds.py b/lib/iris/tests/unit/analysis/cartography/test_rotate_winds.py
new file mode 100644
index 0000000000..ae109e896a
--- /dev/null
+++ b/lib/iris/tests/unit/analysis/cartography/test_rotate_winds.py
@@ -0,0 +1,453 @@
+# (C) British Crown Copyright 2015, 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 the function
+:func:`iris.analysis.cartography.rotate_winds`.
+
+"""
+from __future__ import (absolute_import, division, print_function)
+
+# Import iris.tests first so that some things can be initialised before
+# importing anything else.
+import iris.tests as tests
+
+import numpy as np
+import numpy.ma as ma
+
+import cartopy.crs as ccrs
+from iris.analysis.cartography import rotate_winds, unrotate_pole
+from iris.cube import Cube
+from iris.coords import DimCoord, AuxCoord
+import iris.coord_systems
+
+
+def uv_cubes(x=None, y=None):
+    """Return u, v cubes with a grid in a rotated pole CRS."""
+    cs = iris.coord_systems.RotatedGeogCS(grid_north_pole_latitude=37.5,
+                                          grid_north_pole_longitude=177.5)
+    if x is None:
+        x = np.linspace(311.9, 391.1, 6)
+    if y is None:
+        y = np.linspace(-23.6, 24.8, 5)
+
+    x2d, y2d = np.meshgrid(x, y)
+    u = 10 * (2 * np.cos(2 * np.deg2rad(x2d) + 3 * np.deg2rad(y2d + 30)) ** 2)
+    v = 20 * np.cos(6 * np.deg2rad(x2d))
+    lon = DimCoord(x, standard_name='grid_longitude', units='degrees',
+                   coord_system=cs)
+    lat = DimCoord(y, standard_name='grid_latitude', units='degrees',
+                   coord_system=cs)
+    u_cube = Cube(u, standard_name='x_wind', units='m/s')
+    v_cube = Cube(v, standard_name='y_wind', units='m/s')
+    for cube in (u_cube, v_cube):
+        cube.add_dim_coord(lat.copy(), 0)
+        cube.add_dim_coord(lon.copy(), 1)
+    return u_cube, v_cube
+
+
+def uv_cubes_3d(ref_cube, n_realization=3):
+    """
+    Return 3d u, v cubes with a grid in a rotated pole CRS taken from
+    the provided 2d cube, by adding a realization dimension
+    coordinate bound to teh zeroth dimension.
+
+    """
+    lat = ref_cube.coord('grid_latitude')
+    lon = ref_cube.coord('grid_longitude')
+    x2d, y2d = np.meshgrid(lon.points, lat.points)
+    u = 10 * (2 * np.cos(2 * np.deg2rad(x2d) + 3 * np.deg2rad(y2d + 30)) ** 2)
+    v = 20 * np.cos(6 * np.deg2rad(x2d))
+    # Multiply slices by factor to give variation over 0th dim.
+    factor = np.arange(1, n_realization + 1).reshape(n_realization, 1, 1)
+    u = factor * u
+    v = factor * v
+    realization = DimCoord(np.arange(n_realization), 'realization')
+    u_cube = Cube(u, standard_name='x_wind', units='m/s')
+    v_cube = Cube(v, standard_name='y_wind', units='m/s')
+    for cube in (u_cube, v_cube):
+        cube.add_dim_coord(realization.copy(), 0)
+        cube.add_dim_coord(lat.copy(), 1)
+        cube.add_dim_coord(lon.copy(), 2)
+    return u_cube, v_cube
+
+
+class TestPrerequisites(tests.IrisTest):
+    def test_different_coord_systems(self):
+        u, v = uv_cubes()
+        v.coord('grid_latitude').coord_system = iris.coord_systems.GeogCS(1)
+        with self.assertRaisesRegexp(
+                ValueError, 'Coordinates differ between u and v cubes'):
+            rotate_winds(u, v, iris.coord_systems.OSGB())
+
+    def test_different_xy_coord_systems(self):
+        u, v = uv_cubes()
+        u.coord('grid_latitude').coord_system = iris.coord_systems.GeogCS(1)
+        v.coord('grid_latitude').coord_system = iris.coord_systems.GeogCS(1)
+        with self.assertRaisesRegexp(
+                ValueError,
+                'Coordinate systems of x and y coordinates differ'):
+            rotate_winds(u, v, iris.coord_systems.OSGB())
+
+    def test_different_shape(self):
+        x = np.linspace(311.9, 391.1, 6)
+        y = np.linspace(-23.6, 24.8, 5)
+        u, _ = uv_cubes(x, y)
+        _, v = uv_cubes(x[:-1], y)
+        with self.assertRaisesRegexp(ValueError, 'same shape'):
+            rotate_winds(u, v, iris.coord_systems.OSGB())
+
+    def test_xy_dimensionality(self):
+        u, v = uv_cubes()
+        # Replace 1d lat with 2d lat.
+        x = u.coord('grid_longitude').points
+        y = u.coord('grid_latitude').points
+        x2d, y2d = np.meshgrid(x, y)
+        lat_2d = AuxCoord(y2d, 'grid_latitude', units='degrees',
+                          coord_system=u.coord('grid_latitude').coord_system)
+        for cube in (u, v):
+            cube.remove_coord('grid_latitude')
+            cube.add_aux_coord(lat_2d.copy(), (0, 1))
+
+        with self.assertRaisesRegexp(
+                ValueError,
+                'x and y coordinates must have the same number of dimensions'):
+            rotate_winds(u, v, iris.coord_systems.OSGB())
+
+    def test_dim_mapping(self):
+        x = np.linspace(311.9, 391.1, 3)
+        y = np.linspace(-23.6, 24.8, 3)
+        u, v = uv_cubes(x, y)
+        v.transpose()
+        with self.assertRaisesRegexp(ValueError, 'Dimension mapping'):
+            rotate_winds(u, v, iris.coord_systems.OSGB())
+
+
+class TestAnalyticComparison(tests.IrisTest):
+    @staticmethod
+    def _unrotate_equation(rotated_lons, rotated_lats,
+                           rotated_us, rotated_vs, pole_lon, pole_lat):
+        # Perform a rotated-pole 'unrotate winds' transformation on arrays of
+        # rotated-lat, rotated-lon, u and v.
+        # This can be defined as an analytic function : cf. UMDP015
+
+        # Work out the rotation angles.
+        lambda_angle = np.radians(pole_lon - 180.0)
+        phi_angle = np.radians(90.0 - pole_lat)
+
+        # Get the locations in true lats+lons.
+        trueLongitude, trueLatitude = unrotate_pole(rotated_lons,
+                                                    rotated_lats,
+                                                    pole_lon,
+                                                    pole_lat)
+
+        # Calculate inter-coordinate rotation coefficients.
+        cos_rot = (np.cos(np.radians(rotated_lons)) *
+                   np.cos(np.radians(trueLongitude) - lambda_angle) +
+                   np.sin(np.radians(rotated_lons)) *
+                   np.sin(np.radians(trueLongitude) - lambda_angle) *
+                   np.cos(phi_angle))
+        sin_rot = -((np.sin(np.radians(trueLongitude) - lambda_angle) *
+                     np.sin(phi_angle))
+                    / np.cos(np.radians(rotated_lats)))
+
+        # Matrix-multiply to rotate the vectors.
+        u_true = rotated_us * cos_rot - rotated_vs * sin_rot
+        v_true = rotated_vs * cos_rot + rotated_us * sin_rot
+
+        return u_true, v_true
+
+    def _check_rotated_to_true(self, u_rot, v_rot, target_cs, **kwds):
+        # Run test calculation (numeric).
+        u_true, v_true = rotate_winds(u_rot, v_rot, target_cs)
+
+        # Perform same calculation via the reference method (equations).
+        cs_rot = u_rot.coord('grid_longitude').coord_system
+        pole_lat = cs_rot.grid_north_pole_latitude
+        pole_lon = cs_rot.grid_north_pole_longitude
+        rotated_lons = u_rot.coord('grid_longitude').points
+        rotated_lats = u_rot.coord('grid_latitude').points
+        rotated_lons_2d, rotated_lats_2d = np.meshgrid(
+            rotated_lons, rotated_lats)
+        rotated_u, rotated_v = u_rot.data, v_rot.data
+        u_ref, v_ref = self._unrotate_equation(rotated_lons_2d,
+                                               rotated_lats_2d,
+                                               rotated_u, rotated_v,
+                                               pole_lon, pole_lat)
+
+        # Check that all the numerical results are within given tolerances.
+        self.assertArrayAllClose(u_true.data, u_ref, **kwds)
+        self.assertArrayAllClose(v_true.data, v_ref, **kwds)
+
+    def test_rotated_to_true__small(self):
+        # Check for a small field with varying data.
+        target_cs = iris.coord_systems.GeogCS(6371229)
+        u_rot, v_rot = uv_cubes()
+        self._check_rotated_to_true(u_rot, v_rot, target_cs,
+                                    rtol=1e-5, atol=0.0005)
+
+    def test_rotated_to_true_global(self):
+        # Check for global fields with various constant wind values
+        # - constant in the rotated pole system, that is.
+        # We expect less accuracy where this gets close to the true poles.
+        target_cs = iris.coord_systems.GeogCS(6371229)
+        u_rot, v_rot = uv_cubes(x=np.arange(0, 360.0, 15),
+                                y=np.arange(-89, 89, 10))
+        for vector in ((1, 0), (0, 1), (1, 1), (-3, -1.5)):
+            u_rot.data[...] = vector[0]
+            v_rot.data[...] = vector[1]
+            self._check_rotated_to_true(u_rot, v_rot, target_cs,
+                                        rtol=5e-4, atol=5e-4,
+                                        err_msg='vector={}'.format(vector))
+
+
+class TestRotatedToOSGB(tests.IrisTest):
+    # Define some coordinate ranges for the uv_cubes 'standard' RotatedPole
+    # system, that exceed the OSGB margins, but not by "too much".
+    _rp_x_min, _rp_x_max = -5.0, 5.0
+    _rp_y_min, _rp_y_max = -5.0, 15.0
+
+    def _uv_cubes_limited_extent(self):
+        # Make test cubes suitable for transforming to OSGB, as the standard
+        # 'uv_cubes' result goes too far outside, leading to errors.
+        x = np.linspace(self._rp_x_min, self._rp_x_max, 6)
+        y = np.linspace(self._rp_y_min, self._rp_y_max, 5)
+        return uv_cubes(x=x, y=y)
+
+    def test_name(self):
+        u, v = self._uv_cubes_limited_extent()
+        u.rename('bob')
+        v.rename('alice')
+        ut, vt = rotate_winds(u, v, iris.coord_systems.OSGB())
+        self.assertEqual(ut.name(), 'transformed_' + u.name())
+        self.assertEqual(vt.name(), 'transformed_' + v.name())
+
+    def test_new_coords(self):
+        u, v = self._uv_cubes_limited_extent()
+        x = u.coord('grid_longitude').points
+        y = u.coord('grid_latitude').points
+        x2d, y2d = np.meshgrid(x, y)
+        src_crs = ccrs.RotatedPole(pole_longitude=177.5, pole_latitude=37.5)
+        tgt_crs = ccrs.OSGB()
+        xyz_tran = tgt_crs.transform_points(src_crs, x2d, y2d)
+
+        ut, vt = rotate_winds(u, v, iris.coord_systems.OSGB())
+
+        points = xyz_tran[..., 0].reshape(x2d.shape)
+        expected_x = AuxCoord(points,
+                              standard_name='projection_x_coordinate',
+                              units='m',
+                              coord_system=iris.coord_systems.OSGB())
+        self.assertEqual(ut.coord('projection_x_coordinate'), expected_x)
+        self.assertEqual(vt.coord('projection_x_coordinate'), expected_x)
+
+        points = xyz_tran[..., 1].reshape(y2d.shape)
+        expected_y = AuxCoord(points,
+                              standard_name='projection_y_coordinate',
+                              units='m',
+                              coord_system=iris.coord_systems.OSGB())
+        self.assertEqual(ut.coord('projection_y_coordinate'), expected_y)
+        self.assertEqual(vt.coord('projection_y_coordinate'), expected_y)
+
+    def test_new_coords_transposed(self):
+        u, v = self._uv_cubes_limited_extent()
+        # Transpose cubes so that cube is in xy order rather than the
+        # typical yx order of meshgrid.
+        u.transpose()
+        v.transpose()
+        x = u.coord('grid_longitude').points
+        y = u.coord('grid_latitude').points
+        x2d, y2d = np.meshgrid(x, y)
+        src_crs = ccrs.RotatedPole(pole_longitude=177.5, pole_latitude=37.5)
+        tgt_crs = ccrs.OSGB()
+        xyz_tran = tgt_crs.transform_points(src_crs, x2d, y2d)
+
+        ut, vt = rotate_winds(u, v, iris.coord_systems.OSGB())
+
+        points = xyz_tran[..., 0].reshape(x2d.shape)
+        expected_x = AuxCoord(points,
+                              standard_name='projection_x_coordinate',
+                              units='m',
+                              coord_system=iris.coord_systems.OSGB())
+        self.assertEqual(ut.coord('projection_x_coordinate'), expected_x)
+        self.assertEqual(vt.coord('projection_x_coordinate'), expected_x)
+
+        points = xyz_tran[..., 1].reshape(y2d.shape)
+        expected_y = AuxCoord(points,
+                              standard_name='projection_y_coordinate',
+                              units='m',
+                              coord_system=iris.coord_systems.OSGB())
+        self.assertEqual(ut.coord('projection_y_coordinate'), expected_y)
+        self.assertEqual(vt.coord('projection_y_coordinate'), expected_y)
+        # Check dim mapping for 2d coords is yx.
+        expected_dims = (u.coord_dims('grid_latitude') +
+                         u.coord_dims('grid_longitude'))
+        self.assertEqual(ut.coord_dims('projection_x_coordinate'),
+                         expected_dims)
+        self.assertEqual(ut.coord_dims('projection_y_coordinate'),
+                         expected_dims)
+        self.assertEqual(vt.coord_dims('projection_x_coordinate'),
+                         expected_dims)
+        self.assertEqual(vt.coord_dims('projection_y_coordinate'),
+                         expected_dims)
+
+    def test_orig_coords(self):
+        u, v = self._uv_cubes_limited_extent()
+        ut, vt = rotate_winds(u, v, iris.coord_systems.OSGB())
+        self.assertEqual(u.coord('grid_latitude'), ut.coord('grid_latitude'))
+        self.assertEqual(v.coord('grid_latitude'), vt.coord('grid_latitude'))
+        self.assertEqual(u.coord('grid_longitude'), ut.coord('grid_longitude'))
+        self.assertEqual(v.coord('grid_longitude'), vt.coord('grid_longitude'))
+
+    def test_magnitude_preservation(self):
+        u, v = self._uv_cubes_limited_extent()
+        ut, vt = rotate_winds(u, v, iris.coord_systems.OSGB())
+        orig_sq_mag = u.data**2 + v.data**2
+        res_sq_mag = ut.data**2 + vt.data**2
+        self.assertArrayAllClose(orig_sq_mag, res_sq_mag, rtol=5e-4)
+
+    def test_data_values(self):
+        u, v = self._uv_cubes_limited_extent()
+        # Slice out 4 points that lie in and outside OSGB extent.
+        u = u[1:3, 3:5]
+        v = v[1:3, 3:5]
+        ut, vt = rotate_winds(u, v, iris.coord_systems.OSGB())
+        # Values precalculated and checked.
+        expected_ut_data = np.array([[0.16285514,  0.35323639],
+                                     [1.82650698,  2.62455840]])
+        expected_vt_data = np.array([[19.88979966,  19.01921346],
+                                     [19.88018847,  19.01424281]])
+        # Compare u and v data values against previously calculated values.
+        self.assertArrayAllClose(ut.data, expected_ut_data, rtol=1e-5)
+        self.assertArrayAllClose(vt.data, expected_vt_data, rtol=1e-5)
+
+    def test_nd_data(self):
+        u2d, y2d = self._uv_cubes_limited_extent()
+        u, v = uv_cubes_3d(u2d)
+        u = u[:, 1:3, 3:5]
+        v = v[:, 1:3, 3:5]
+        ut, vt = rotate_winds(u, v, iris.coord_systems.OSGB())
+        # Values precalculated and checked (as test_data_values above),
+        # then scaled by factor [1, 2, 3] along 0th dim (see uv_cubes_3d()).
+        expected_ut_data = np.array([[0.16285514,  0.35323639],
+                                     [1.82650698,  2.62455840]])
+        expected_vt_data = np.array([[19.88979966,  19.01921346],
+                                     [19.88018847,  19.01424281]])
+        factor = np.array([1, 2, 3]).reshape(3, 1, 1)
+        expected_ut_data = factor * expected_ut_data
+        expected_vt_data = factor * expected_vt_data
+        # Compare u and v data values against previously calculated values.
+        self.assertArrayAlmostEqual(ut.data, expected_ut_data)
+        self.assertArrayAlmostEqual(vt.data, expected_vt_data)
+
+    def test_transposed(self):
+        # Test case where the coordinates are not ordered yx in the cube.
+        u, v = self._uv_cubes_limited_extent()
+        # Slice out 4 points that lie in and outside OSGB extent.
+        u = u[1:3, 3:5]
+        v = v[1:3, 3:5]
+        # Transpose cubes (in-place)
+        u.transpose()
+        v.transpose()
+        ut, vt = rotate_winds(u, v, iris.coord_systems.OSGB())
+        # Values precalculated and checked.
+        expected_ut_data = np.array([[0.16285514,  0.35323639],
+                                     [1.82650698,  2.62455840]]).T
+        expected_vt_data = np.array([[19.88979966,  19.01921346],
+                                     [19.88018847,  19.01424281]]).T
+        # Compare u and v data values against previously calculated values.
+        self.assertArrayAllClose(ut.data, expected_ut_data, rtol=1e-5)
+        self.assertArrayAllClose(vt.data, expected_vt_data, rtol=1e-5)
+
+
+class TestMasking(tests.IrisTest):
+    def test_rotated_to_osgb(self):
+        # Rotated Pole data with large extent.
+        x = np.linspace(311.9, 391.1, 10)
+        y = np.linspace(-23.6, 24.8, 8)
+        u, v = uv_cubes(x, y)
+        ut, vt = rotate_winds(u, v, iris.coord_systems.OSGB())
+
+        # Ensure cells with discrepancies in magnitude are masked.
+        self.assertTrue(ma.isMaskedArray(ut.data))
+        self.assertTrue(ma.isMaskedArray(vt.data))
+
+        # Snapshot of mask with fixed tolerance of atol=2e-3
+        expected_mask = np.array([[1, 1, 1, 0, 0, 0, 0, 0, 0, 1],
+                                  [1, 1, 1, 0, 0, 0, 0, 0, 0, 1],
+                                  [1, 1, 1, 1, 0, 0, 0, 0, 1, 1],
+                                  [1, 1, 1, 1, 0, 0, 0, 0, 1, 1],
+                                  [1, 1, 1, 1, 0, 0, 0, 0, 1, 1],
+                                  [1, 1, 1, 1, 1, 0, 0, 1, 1, 1],
+                                  [1, 1, 1, 1, 1, 0, 0, 1, 1, 1],
+                                  [1, 1, 1, 1, 1, 0, 0, 1, 1, 1]], np.bool)
+        self.assertArrayEqual(expected_mask, ut.data.mask)
+        self.assertArrayEqual(expected_mask, vt.data.mask)
+
+        # Check unmasked values have sufficiently small error in mag.
+        expected_mag = np.sqrt(u.data**2 + v.data**2)
+        # Use underlying data to ignore mask in calculation.
+        res_mag = np.sqrt(ut.data.data**2 + vt.data.data**2)
+        # Calculate percentage error (note there are no zero magnitudes
+        # so we can divide safely).
+        anom = 100.0 * np.abs(res_mag - expected_mag) / expected_mag
+        self.assertTrue(anom[~ut.data.mask].max() < 0.1)
+
+    def test_rotated_to_unrotated(self):
+        # Suffiently accurate so that no mask is introduced.
+        u, v = uv_cubes()
+        ut, vt = rotate_winds(u, v, iris.coord_systems.GeogCS(6371229))
+        self.assertFalse(ma.isMaskedArray(ut.data))
+        self.assertFalse(ma.isMaskedArray(vt.data))
+
+
+class TestRoundTrip(tests.IrisTest):
+    def test_rotated_to_unrotated(self):
+        # Check ability to use 2d coords as input.
+        u, v = uv_cubes()
+        ut, vt = rotate_winds(u, v, iris.coord_systems.GeogCS(6371229))
+        # Remove  grid lat and lon, leaving 2d projection coords.
+        ut.remove_coord('grid_latitude')
+        vt.remove_coord('grid_latitude')
+        ut.remove_coord('grid_longitude')
+        vt.remove_coord('grid_longitude')
+        # Change back.
+        orig_cs = u.coord('grid_latitude').coord_system
+        res_u, res_v = rotate_winds(ut, vt, orig_cs)
+        # Check data values - limited accuracy due to numerical approx.
+        self.assertArrayAlmostEqual(res_u.data, u.data, decimal=3)
+        self.assertArrayAlmostEqual(res_v.data, v.data, decimal=3)
+        # Check coords locations.
+        x2d, y2d = np.meshgrid(u.coord('grid_longitude').points,
+                               u.coord('grid_latitude').points)
+        # Shift longitude from 0 to 360 -> -180 to 180.
+        x2d = np.where(x2d > 180, x2d - 360, x2d)
+        res_x = res_u.coord('projection_x_coordinate',
+                            coord_system=orig_cs).points
+        res_y = res_u.coord('projection_y_coordinate',
+                            coord_system=orig_cs).points
+        self.assertArrayAlmostEqual(res_x, x2d)
+        self.assertArrayAlmostEqual(res_y, y2d)
+        res_x = res_v.coord('projection_x_coordinate',
+                            coord_system=orig_cs).points
+        res_y = res_v.coord('projection_y_coordinate',
+                            coord_system=orig_cs).points
+        self.assertArrayAlmostEqual(res_x, x2d)
+        self.assertArrayAlmostEqual(res_y, y2d)
+
+
+if __name__ == "__main__":
+    tests.main()