diff --git a/CHANGES b/CHANGES index 879b3c2dc5..a266221cb2 100644 --- a/CHANGES +++ b/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 -------------------- diff --git a/docs/iris/src/whats_new.rst b/docs/iris/src/whats_new.rst index 8260dc12d2..07a2d6b218 100644 --- a/docs/iris/src/whats_new.rst +++ b/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 diff --git a/lib/iris/analysis/cartography.py b/lib/iris/analysis/cartography.py index 6c18fd8062..8898d11347 100644 --- a/lib/iris/analysis/cartography.py +++ b/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 diff --git a/lib/iris/tests/results/analysis/project/default_source_cs.cml b/lib/iris/tests/results/analysis/project/default_source_cs.cml new file mode 100644 index 0000000000..998a174278 --- /dev/null +++ b/lib/iris/tests/results/analysis/project/default_source_cs.cml @@ -0,0 +1,104 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/lib/iris/tests/results/visual_tests/test_analysis.TestProject.test_cartopy_projection.0.png b/lib/iris/tests/results/visual_tests/test_analysis.TestProject.test_cartopy_projection.0.png new file mode 100644 index 0000000000..13c2e3e18e Binary files /dev/null and b/lib/iris/tests/results/visual_tests/test_analysis.TestProject.test_cartopy_projection.0.png differ diff --git a/lib/iris/tests/test_analysis.py b/lib/iris/tests/test_analysis.py index a15d8939dc..58ad856fb7 100644 --- a/lib/iris/tests/test_analysis.py +++ b/lib/iris/tests/test_analysis.py @@ -20,10 +20,13 @@ # import iris tests first so that some things can be initialised before importing anything else import iris.tests as tests +import itertools import os import unittest import zlib +import cartopy.crs as ccrs +import matplotlib import matplotlib.pyplot as plt import numpy import shapely.geometry @@ -558,6 +561,89 @@ def test_shared_xy(self): self.assertTrue(numpy.allclose(weights, target)) +class TestProject(tests.GraphicsTest): + def setUp(self): + cube = iris.tests.stock.realistic_4d_no_derived() + # Remove some slices to speed testing. + self.cube = cube[0:2, 0:3] + self.target_proj = ccrs.Robinson() + + def test_bad_resolution(self): + with self.assertRaises(ValueError): + iris.analysis.cartography.project(self.cube, + self.target_proj, + nx=-200, ny=200) + with self.assertRaises(ValueError): + iris.analysis.cartography.project(self.cube, + self.target_proj, + nx=200, ny='abc') + + def test_missing_latlon(self): + cube = self.cube.copy() + cube.remove_coord('grid_latitude') + with self.assertRaises(ValueError): + iris.analysis.cartography.project(cube, self.target_proj) + cube = self.cube.copy() + cube.remove_coord('grid_longitude') + with self.assertRaises(ValueError): + iris.analysis.cartography.project(cube, self.target_proj) + self.cube.remove_coord('grid_longitude') + self.cube.remove_coord('grid_latitude') + with self.assertRaises(ValueError): + iris.analysis.cartography.project(self.cube, self.target_proj) + + def test_default_resolution(self): + new_cube, extent = iris.analysis.cartography.project(self.cube, + self.target_proj) + self.assertEqual(new_cube.shape, self.cube.shape) + + @iris.tests.skip_data + def test_cartopy_projection(self): + cube = iris.load_strict(tests.get_data_path(('PP', 'aPPglob1', + 'global.pp'))) + projections = {} + projections['RotatedPole'] = ccrs.RotatedPole(pole_longitude=177.5, + pole_latitude=37.5) + projections['Robinson'] = ccrs.Robinson() + projections['PlateCarree'] = ccrs.PlateCarree() + projections['NorthPolarStereo'] = ccrs.NorthPolarStereo() + projections['Orthographic'] = ccrs.Orthographic(central_longitude=-90, + central_latitude=45) + projections['InterruptedGoodeHomolosine'] = ccrs.InterruptedGoodeHomolosine() + projections['LambertCylindrical'] = ccrs.LambertCylindrical() + + # Set up figure + fig = plt.figure(figsize=(10, 10)) + gs = matplotlib.gridspec.GridSpec(nrows=3, ncols=3, hspace=1.5, wspace=0.5) + for subplot_spec, (name, target_proj) in itertools.izip(gs, projections.iteritems()): + # Set up axes and title + ax = plt.subplot(subplot_spec, frameon=False, projection=target_proj) + ax.set_title(name) + # Transform cube to target projection + new_cube, extent = iris.analysis.cartography.project(cube, target_proj, + nx=150, ny=150) + # Plot + plt.pcolor(new_cube.coord('projection_x_coordinate').points, + new_cube.coord('projection_y_coordinate').points, + new_cube.data) + # Add coastlines + ax.coastlines() + + # Tighten up layout + gs.tight_layout(plt.gcf()) + + # Verify resulting plot + self.check_graphic() + + @iris.tests.skip_data + def test_no_coord_system(self): + cube = iris.load_strict(tests.get_data_path(('PP', 'aPPglob1', 'global.pp'))) + cube.coord('longitude').coord_system = None + cube.coord('latitude').coord_system = None + new_cube, extent = iris.analysis.cartography.project(cube, self.target_proj) + self.assertCML(new_cube, ('analysis', 'project', 'default_source_cs.cml')) + + if __name__ == "__main__": - unittest.main() - + tests.main() +