Added cube project function to iris.analysis.cartography

CHANGES

@@ -14,6 +14,8 @@ Features added
 * A more explicit set of load functions, which also allow the automatic
   cube merging to be bypassed as a last resort.
 * Save netCDF files with an unlimited dimension.
+* The ability to project a cube with a lat-lon or rotated lat-lon coordinate
+  system into a range of map projections e.g. Polar Stereographic.
 
 Incompatible changes
 --------------------

docs/iris/src/whats_new.rst

@@ -41,6 +41,9 @@ A summary of the main features added with version 1.0:
 * Save netCDF files with an unlimited dimension.
 * A more explicit set of load functions, which also allow the automatic
   cube merging to be bypassed as a last resort.
+* The ability to project a cube with a lat-lon or rotated lat-lon coordinate
+  system into a range of map projections e.g. Polar Stereographic.
+
 
 Incompatible changes
 --------------------
@@ -229,6 +232,44 @@ now use the :func:`iris.load_cube()` and :func:`iris.load_cubes()`
 functions instead.
 
 
+Cube projection
+===============
+
+Iris now has the ability to project a cube into a number of map projections.
+This functionality is provided by :func:`iris.analysis.cartography.project()`.
+For example::
+
+    import iris
+    import cartopy
+    import matplotlib.pyplot as plt
+
+    # Load data
+    cube = iris.load_cube(iris.sample_data_path('air_temp.pp'))
+
+    # Transform cube to target projection
+    target_proj = cartopy.crs.RotatedPole(pole_longitude=177.5,
+                                          pole_latitude=37.5)
+    new_cube, extent = iris.analysis.cartography.project(cube, target_proj)
+
+    # Plot
+    plt.axes(projection=target_proj)
+    plt.pcolor(new_cube.coord('projection_x_coordinate').points,
+               new_cube.coord('projection_y_coordinate').points,
+               new_cube.data)
+    plt.gca().coastlines()
+    plt.show()
+
+This function is intended to be used in cases where the cube's coordinates
+prevent one from directly visualising the data, e.g. when the longitude
+and latitude are two dimensional and do not make up a regular grid. The
+function uses a nearest neighbour approach rather than any form of
+linear/non-linear interpolation to determine the data value of each cell
+in the resulting cube. Consequently it may have an adverse effect on the
+statistics of the data e.g. the mean and standard deviation will not be
+preserved. This function currently assumes global data and will if
+necessary extrapolate beyond the geographical extent of the source cube.
+
+
 Other changes
 =============
 * Cube summaries are now more readable when the scalar coordinates

lib/iris/analysis/cartography.py

@@ -18,12 +18,12 @@
 Various utilities and numeric transformations relevant to cartography.
 
 """
-
+import itertools
 import math
 import warnings
 
+import cartopy.img_transform
 import numpy
-import cartopy.crs
 
 import iris.analysis
 import iris.coords
@@ -332,3 +332,230 @@ def area_weights(cube):
     broad_weights = ll_weights * other_dims_array
 
     return broad_weights
+
+
+def project(cube, target_proj, nx=None, ny=None):
+    """
+    Return a new cube that is the result of projecting a cube from its
+    coordinate system into a specified projection e.g. Robinson or Polar
+    Stereographic. This function is intended to be used in cases where the
+    cube's coordinates prevent one from directly visualising the data, e.g.
+    when the longitude and latitude are two dimensional and do not make up
+    a regular grid.
+
+    Args:
+        * cube
+            An instance of :class:`iris.cube.Cube`.
+        * target_proj
+            An instance of the Cartopy Projection class, or an instance of
+            :class:`iris.coord_systems.CoordSystem` from which a projection
+            will be obtained.
+    Kwargs:
+        * nx
+            Desired number of sample points in the x direction.
+        * ny
+            Desired number of sample points in the y direction.
+
+    Returns:
+        An instance of :class:`iris.cube.Cube` and a list describing the
+        extent of the projection.
+
+    .. note::
+        This function uses a nearest neighbour approach rather than any form
+        of linear/non-linear interpolation to determine the data value of each
+        cell in the resulting cube. Consequently it may have an adverse effect
+        on the statistics of the data e.g. the mean and standard deviation
+        will not be preserved.
+
+    .. note::
+        This function assumes global data and will if necessary extrapolate
+        beyond the geographical extent of the source cube using a nearest
+        neighbour approach.
+
+    """
+    try:
+        lat_coord, lon_coord = _get_lat_lon_coords(cube)
+    except IndexError:
+        raise ValueError('Cannot get latitude/longitude ' \
+                         'coordinates from cube {!r}.'.format(cube.name()))
+
+    if lat_coord.coord_system != lon_coord.coord_system:
+        raise ValueError('latitude and longitude coords appear to have' \
+                         'different coordinates systems.')
+
+    if lon_coord.units != 'degrees':
+        lon_coord = lon_coord.unit_converted('degrees')
+    if lat_coord.units != 'degrees':
+        lat_coord = lat_coord.unit_converted('degrees')
+
+    # Determine source coordinate system
+    if lat_coord.coord_system is None:
+        # Assume WGS84 latlon if unspecified
+        warnings.warn('Coordinate system of latitude and longitude '\
+                      'coordinates is not specified. Assuming WGS84 Geodetic.')
+        orig_cs = iris.coord_systems.GeogCS(semi_major_axis=6378137.0,
+                                            inverse_flattening=298.257223563)
+    else:
+        orig_cs = lat_coord.coord_system
+
+    # Convert to cartopy crs
+    source_cs = orig_cs.as_cartopy_crs()
+
+    # Obtain coordinate arrays (ignoring bounds) and convert to 2d
+    # if not already.
+    source_x = lon_coord.points
+    source_y = lat_coord.points
+    if source_x.ndim != 2 or source_y.ndim != 2:
+        source_x, source_y = numpy.meshgrid(source_x, source_y)
+
+    # Calculate target grid
+    if isinstance(target_proj, iris.coord_systems.CoordSystem):
+        target_proj = target_proj.as_cartopy_projection()
+
+    # Resolution of new grid
+    if nx == None:
+        nx = source_x.shape[1]
+    if ny == None:
+        ny = source_x.shape[0]
+
+    target_x, target_y, extent = cartopy.img_transform.mesh_projection(target_proj,
+                                                                       nx, ny)
+
+    # Determine dimension mappings - expect either 1d or 2d
+    if lat_coord.ndim != lon_coord.ndim:
+        raise ValueError("The latitude and longitude coordinates have "
+                         "different dimensionality.")
+
+    latlon_ndim = lat_coord.ndim
+    lon_dims = cube.coord_dims(lon_coord)
+    lat_dims = cube.coord_dims(lat_coord)
+
+    if latlon_ndim == 1:
+        xdim = lon_dims[0]
+        ydim = lat_dims[0]
+    elif latlon_ndim == 2:
+        if lon_dims != lat_dims:
+            raise ValueError("The 2d latitude and longitude coordinates "
+                             "correspond to different dimensions.")
+        # If coords are 2d assume that grid is ordered such that x corresponds
+        # to the last dimension (shortest stride).
+        xdim = lon_dims[1]
+        ydim = lon_dims[0]
+    else:
+        raise ValueError('Expected the latitude and longitude coordinates '\
+                         'to have 1 or 2 dimensions, got {} and '\
+                         '{}.'.format(lat_coord.ndim, lon_coord.ndim))
+
+    # Create array to store regridded data
+    new_shape = list(cube.shape)
+    new_shape[xdim] = nx
+    new_shape[ydim] = ny
+    new_data = numpy.ma.zeros(new_shape, cube.data.dtype)
+
+    # Create iterators to step through cube data in lat long slices
+    new_shape[xdim] = 1
+    new_shape[ydim] = 1
+    index_it = numpy.ndindex(*new_shape)
+    if lat_coord.ndim == 1 and lon_coord.ndim == 1:
+        slice_it = cube.slices([lat_coord, lon_coord])
+    elif lat_coord.ndim == 2 and lon_coord.ndim == 2:
+        slice_it = cube.slices(lat_coord)
+    else:
+        raise ValueError('Expected the latitude and longitude coordinates '\
+                         'to have 1 or 2 dimensions, got {} and '\
+                         '{}.'.format(lat_coord.ndim, lon_coord.ndim))
+
+    ## Mask out points outside of extent in source_cs - disabled until
+    ## a way to specify global/limited extent is agreed upon and code
+    ## is generalised to handle -180 to +180, 0 to 360 and >360 longitudes.
+    #source_desired_xy = source_cs.transform_points(target_proj,
+    #                                               target_x.flatten(),
+    #                                               target_y.flatten())
+    #if numpy.any(source_x < 0.0) and numpy.any(source_x > 180.0):
+    #    raise ValueError('Unable to handle range of longitude.')
+    ## This does not work in all cases e.g. lon > 360
+    #if numpy.any(source_x > 180.0):
+    #    source_desired_x = (source_desired_xy[:, 0].reshape(ny, nx) + 360.0) % 360.0
+    #else:
+    #    source_desired_x = source_desired_xy[:, 0].reshape(ny, nx)
+    #source_desired_y = source_desired_xy[:, 1].reshape(ny, nx)
+    #outof_extent_points = ((source_desired_x < source_x.min()) |
+    #                       (source_desired_x > source_x.max()) |
+    #                       (source_desired_y < source_y.min()) |
+    #                       (source_desired_y > source_y.max()))
+    ## Make array a mask by default (rather than a single bool) to allow mask to be
+    ## assigned to slices.
+    #new_data.mask = numpy.zeros(new_shape)
+
+    # Step through cube data, regrid onto desired projection and insert results
+    # in new_data array
+    for index, ll_slice in itertools.izip(index_it, slice_it):
+        # Regrid source data onto target grid
+        index = list(index)
+        index[xdim] = slice(None, None)
+        index[ydim] = slice(None, None)
+        new_data[index] = cartopy.img_transform.regrid(ll_slice.data,
+                                                       source_x, source_y,
+                                                       source_cs, target_proj,
+                                                       target_x, target_y)
+
+        ## Mask out points beyond extent
+        #new_data[index].mask[outof_extent_points] = True
+
+    # Remove mask if it is unnecessary
+    if not numpy.any(new_data.mask):
+        new_data = new_data.data
+
+    # Create new cube
+    new_cube = iris.cube.Cube(new_data)
+
+    # Add new grid coords
+    x_coord = iris.coords.DimCoord(target_x[0, :], 'projection_x_coordinate')
+    y_coord = iris.coords.DimCoord(target_y[:, 0], 'projection_y_coordinate')
+    new_cube.add_dim_coord(x_coord, xdim)
+    new_cube.add_dim_coord(y_coord, ydim)
+
+    # Add resampled lat/lon in original coord system
+    source_desired_xy = source_cs.transform_points(target_proj,
+                                                   target_x.flatten(),
+                                                   target_y.flatten())
+    new_lon_points = source_desired_xy[:, 0].reshape(ny, nx)
+    new_lat_points = source_desired_xy[:, 1].reshape(ny, nx)
+    new_lon_coord = iris.coords.AuxCoord(new_lon_points,
+                                         standard_name='longitude',
+                                         units='degrees',
+                                         coord_system=orig_cs)
+    new_lat_coord = iris.coords.AuxCoord(new_lat_points,
+                                         standard_name='latitude',
+                                         units='degrees',
+                                         coord_system=orig_cs)
+    new_cube.add_aux_coord(new_lon_coord, [ydim, xdim])
+    new_cube.add_aux_coord(new_lat_coord, [ydim, xdim])
+
+    coords_to_ignore = set()
+    coords_to_ignore.update(cube.coords(contains_dimension=xdim))
+    coords_to_ignore.update(cube.coords(contains_dimension=ydim))
+    for coord in cube.dim_coords:
+        if coord not in coords_to_ignore:
+            new_cube.add_dim_coord(coord.copy(), cube.coord_dims(coord))
+    for coord in cube.aux_coords:
+        if coord not in coords_to_ignore:
+            new_cube.add_aux_coord(coord.copy(), cube.coord_dims(coord))
+    discarded_coords = coords_to_ignore.difference([lat_coord, lon_coord])
+    if discarded_coords:
+        warnings.warn('Discarding coordinates that share dimensions with ' \
+                      '{} and {}: {}'.format(lat_coord.name(),
+                                             lon_coord.name(),
+                                             [coord.name() for
+                                              coord in discarded_coords]))
+
+    # TODO handle derived coords/aux_factories
+
+    # Copy metadata across
+    new_cube.metadata = cube.metadata
+
+    # Record transform in cube's history
+    new_cube.add_history('Converted from {} to {}'.format(type(source_cs).__name__,
+                                                          type(target_proj).__name__))
+
+    return new_cube, extent