-
Notifications
You must be signed in to change notification settings - Fork 300
ENH: Added coordinate inversion functionality #1731
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Changes from all commits
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -35,6 +35,53 @@ | |
| from iris.coord_systems import GeogCS, RotatedGeogCS | ||
|
|
||
|
|
||
| def geodetic(shape, with_bounds=True, xlim=None, ylim=None, circular=False): | ||
| """ | ||
| Return a global lat-lon cube of specified shape. | ||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This could do with adding "2D", clarifying "global", and describing |
||
|
|
||
| """ | ||
| if xlim is None: | ||
| xlim = (-180, 180) | ||
| if ylim is None: | ||
| ylim = (-90, 90) | ||
| shape = np.array(shape) | ||
| assert shape.ndim == 1 | ||
| assert shape.shape[0] == 2 | ||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I find the conversion to an array and subsequent assertions somewhat confusing ... how about just: |
||
|
|
||
| data = np.arange(np.product(shape)).reshape(shape) | ||
| data = data.astype('int32') | ||
|
|
||
| if xlim[0] > xlim[1]: | ||
| data = data[:, ::-1] | ||
| if ylim[0] > ylim[1]: | ||
| data = data[::-1] | ||
|
|
||
| cube = iris.cube.Cube(data) | ||
| crs = iris.coord_systems.GeogCS(6371229) | ||
|
|
||
| x_bound = np.linspace(xlim[0], xlim[1], endpoint=True, num=shape[1]+1) | ||
| x_bounds = np.array([x_bound[:-1], x_bound[1:]]).T | ||
| x_points = x_bounds.mean(axis=1) | ||
|
|
||
| y_bound = np.linspace(ylim[0], ylim[1], endpoint=True, num=shape[0]+1) | ||
| y_bounds = np.array([y_bound[:-1], y_bound[1:]]).T | ||
| y_points = y_bounds.mean(axis=1) | ||
|
|
||
| if not with_bounds: | ||
| x_bounds = None | ||
| y_bounds = None | ||
| x_coord = iris.coords.DimCoord(x_points, standard_name='longitude', | ||
| units='degree_east', | ||
| coord_system=crs, bounds=x_bounds, | ||
| circular=circular) | ||
| y_coord = iris.coords.DimCoord(y_points, standard_name='latitude', | ||
| units='degree_north', | ||
| coord_system=crs, bounds=y_bounds) | ||
| cube.add_dim_coord(y_coord, 0) | ||
| cube.add_dim_coord(x_coord, 1) | ||
| return cube | ||
|
|
||
|
|
||
| def lat_lon_cube(): | ||
| """ | ||
| Returns a cube with a latitude and longitude suitable for testing | ||
|
|
||
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,68 @@ | ||
| # (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/>. | ||
| """Test function :func:`iris.util.invert_coordinate`.""" | ||
|
|
||
| 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 | ||
|
|
||
| import unittest | ||
|
|
||
| import iris | ||
| import iris.tests.stock as stock | ||
|
|
||
|
|
||
| class TestAll(iris.tests.IrisTest): | ||
| def setUp(self): | ||
| self.cube = stock.geodetic((2, 2)) | ||
|
|
||
| def test_inversion_values(self): | ||
| tar = self.cube.data[::-1] | ||
|
|
||
| iris.util.invert_coordinate(self.cube, 'latitude') | ||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Really you should explicitly import |
||
|
|
||
| self.assertArrayAlmostEqual(self.cube.data, tar) | ||
|
|
||
| tar = [45., -45.] | ||
| self.assertArrayAlmostEqual(self.cube.coord('latitude').points, tar) | ||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Please split this out as a separate test method. |
||
|
|
||
| tar = [[90, 0], [0, -90]] | ||
| self.assertArrayAlmostEqual(self.cube.coord('latitude').bounds, tar) | ||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Please split this out as a separate test method. |
||
|
|
||
| def test_aux_coord(self): | ||
| coord = iris.coords.AuxCoord.from_coord(self.cube.coord('latitude')) | ||
| self.cube.remove_coord('latitude') | ||
| self.cube.add_aux_coord(coord, 0) | ||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. There an |
||
| msg = ('Only an inversion of a dimension coordinate is supported ' | ||
| '\(latitude\)') | ||
| with self.assertRaisesRegexp(RuntimeError, msg): | ||
| iris.util.invert_coordinate(self.cube, 'latitude') | ||
|
|
||
| def test_double_inversion(self): | ||
| # Check that only the intended changes are made by double inversion. | ||
| # We should have exactly the same thing as we started with. | ||
| tar = self.cube.copy() | ||
| iris.util.invert_coordinate(self.cube, 'latitude') | ||
| iris.util.invert_coordinate(self.cube, 'latitude') | ||
| self.assertEqual(self.cube, tar) | ||
|
|
||
|
|
||
| if __name__ == '__main__': | ||
| unittest.main() | ||
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -41,6 +41,54 @@ | |
| import iris.exceptions | ||
|
|
||
|
|
||
| def _cube_cubes_arg(cubes): | ||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This doesn't really do enough to warrant existing as a separate function. |
||
| """ | ||
| Function for ensuring that we return a Cubelist, irrespective of whether | ||
| a Cube or a CubeList has been provided. | ||
|
|
||
| Args: | ||
|
|
||
| cubes (:class:`iris.cube.Cube` or :class:`iris.cube.CubeList`) | ||
|
|
||
| Returns: | ||
|
|
||
| :class:`iris.cube.CubeList` | ||
|
|
||
| """ | ||
| # Function for ensuring that we always deal with a cubelists | ||
| if isinstance(cubes, iris.cube.Cube): | ||
| cubes = iris.cube.CubeList([cubes]) | ||
| return cubes | ||
|
|
||
|
|
||
| def invert_coordinate(cubes, coordinate): | ||
| """ | ||
| In-place operation to reverse the direction of a coordinate in a cube or | ||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. If this is going to be limited to dimension coordinates then it should say so. Have you considered allowing it to operate on 1-dimensional auxiliary coordinates as well?
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Re. "In-place" ... do you need to be in-place for memory reasons? If not, you could make it more robust (and address @ajdawson's comment) by just indexing the cube, e.g. |
||
| list of cubes. | ||
|
|
||
| Args: | ||
|
|
||
| * cubes (:class:`iris.cube.Cube` or :class:`iris.cube.CubeList`): | ||
| Cube(s) to have their coordinate inverted | ||
|
|
||
| * coordinate ('string' or :class:`iris.coord.Coord`): | ||
|
|
||
| """ | ||
| cubes = _cube_cubes_arg(cubes) | ||
| for cube in cubes: | ||
| coord = cube.coord(coordinate) | ||
| if coord not in cube.dim_coords: | ||
| raise RuntimeError('Only an inversion of a dimension coordinate ' | ||
| 'is supported ({})'.format(coord.name())) | ||
| dims = cube.coord_dims(coord)[0] | ||
|
|
||
| cube.data = reverse(cube.data, dims) | ||
| coord.points = reverse(coord.points, 0) | ||
|
|
||
| if coord.has_bounds(): | ||
| coord.bounds = reverse(coord.bounds, [0, 1]) | ||
|
|
||
|
|
||
| def broadcast_weights(weights, array, dims): | ||
| """ | ||
| Broadcast a weights array to the shape of another array. | ||
|
|
||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There's already a
lat_lon_cubefunction in this module which does a very similar (albeit more limited) job. It'd be nicer if we could just extendlat_lon_cube.