Skip to content
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.

Sign up for GitHub

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

Make difference module handle circular cubes #1990

Closed
wants to merge 24 commits into from
Closed

Make difference module handle circular cubes #1990

Show file tree
Hide file tree
Changes from 16 commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
50d08bb
Resolved merge conflict
MO-PeterJordan Apr 22, 2024
e4bce85
Added another failing test for difference module. Differences calcula…
MO-PeterJordan Apr 22, 2024
beabc9e
Resolved merge conflicts
MO-PeterJordan Apr 23, 2024
141ee1f
Auto-formatting
MO-PeterJordan Apr 23, 2024
9535d45
Import Bugfix
MO-PeterJordan Apr 23, 2024
c442023
Ran isort against spatial.py
MO-PeterJordan Apr 23, 2024
62dc4d9
Added test for appropriate error for unhandled cubes and made code ma…
MO-PeterJordan Apr 23, 2024
636ec69
iSort on spatial.py again
MO-PeterJordan Apr 23, 2024
fa2b8b0
Another attempt at autoformatting
MO-PeterJordan Apr 23, 2024
53659a6
Tidied up
MO-PeterJordan Apr 23, 2024
617bd27
Removed confusing variable name
MO-PeterJordan Apr 24, 2024
0910cd6
Autoformatting
MO-PeterJordan Apr 24, 2024
8082dd0
Tried manually sorting import order as isort still failing in github …
MO-PeterJordan Apr 24, 2024
22e920a
And again
MO-PeterJordan Apr 24, 2024
ce71d4c
And again
MO-PeterJordan Apr 24, 2024
b74833b
Linting suggestion
MO-PeterJordan Apr 24, 2024
ccdc67f
Review Actions
MO-PeterJordan May 20, 2024
c51ad05
Ran isort to fix gha
MO-PeterJordan May 23, 2024
8eb9195
Added test for cubes with flipped axes (lonitude increasing along col…
MO-PeterJordan Jun 17, 2024
939ac46
Merge branch 'master' into pjordan-make-difference-module-handle-circ…
MO-PeterJordan Jun 17, 2024
e3496eb
black, isort, flake8
MO-PeterJordan Jun 17, 2024
175e825
Ran isort in isolation
MO-PeterJordan Jun 17, 2024
6fbc3eb
Fixed test failures resulting from breaking changes made by upstream PR
MO-PeterJordan Jun 17, 2024
31fceb3
Black
MO-PeterJordan Jun 17, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 30 additions & 11 deletions improver/utilities/spatial.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,8 @@
import numpy as np
from cartopy.crs import CRS
from cf_units import Unit
from iris.coords import AuxCoord, CellMethod
from iris.coord_systems import GeogCS
from iris.coords import AuxCoord, CellMethod, Coord
from iris.cube import Cube, CubeList
from numpy import ndarray
from numpy.ma import MaskedArray
Expand Down Expand Up @@ -190,6 +191,10 @@ class DifferenceBetweenAdjacentGridSquares(BasePlugin):
individually.
"""

@staticmethod
def _axis_wraps_around_meridian(axis: Coord, cube: Cube):
return axis.circular and axis == cube.coord(axis="x")

@staticmethod
def _update_metadata(diff_cube: Cube, coord_name: str, cube_name: str) -> None:
"""Rename cube, add attribute and cell method to describe difference.
Expand All @@ -209,9 +214,8 @@ def _update_metadata(diff_cube: Cube, coord_name: str, cube_name: str) -> None:
diff_cube.attributes["form_of_difference"] = "forward_difference"
diff_cube.rename("difference_of_" + cube_name)

@staticmethod
def create_difference_cube(
cube: Cube, coord_name: str, diff_along_axis: ndarray
self, cube: Cube, coord_name: str, diff_along_axis: ndarray
) -> Cube:
"""
Put the difference array into a cube with the appropriate
Expand All @@ -230,8 +234,17 @@ def create_difference_cube(
Cube containing the differences calculated along the
specified axis.
"""
points = cube.coord(coord_name).points
axis = cube.coord(coord_name)
points = axis.points
mean_points = (points[1:] + points[:-1]) / 2
if self._axis_wraps_around_meridian(axis, cube):
if type(axis.coord_system) != GeogCS:
raise ValueError(
"DifferenceBetweenAdjacentGridSquares does not currently support cubes with "
"circular x-axis that do not use a geographic (i.e. latlon) coordinate system."
)
extra_mean_point = np.mean([points[-1], (points[0] + 360)]) % 360
mean_points = np.hstack([mean_points, extra_mean_point])
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Question for the reviewer:

This will only work if the 'wrap around' point is positive. This is the case for test cubes that are 'symmetric' in nature (eg. points in cube at [120, 0, 120] or [-175, 165, 155, ..., -25, -15, -5, 5, 15, 25, ..., 165, 175] - wrap around point = 180) but not ones that are 'offset' (eg. [-170, -160, -150, ..., -20, -10, 0, 10, 20, ..., 150, 160, 170, 180] - wrap around point = -175).

In the offset case given above, this would calculate the 'wrap around point' to be 185. Not 'wrong' per-se but certainly potentially confusing.

I could add extra logic to make this work for offset grids. However given that global data from STaGE is 'symmetrical', and the cubes produced by the set_up_variable_cube test helper match this convention, I'm not sure whether or not this is something worth investing time in.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

...Thinking about it, for the symmetrical case above, the np.mean(...) line always resolves to 180, so I could get rid of it altogether and just always put the wrap around point at 180.

The above does better handle the case where the cube uses 0 -> 360 instead of -180 -> 180 but maybe this never comes up.

I suppose the line could also be used as a test. I.e. if np.mean(...) != 180 raise ValueError("Difference Plugin only handles circular cubes if their points are spread symmetrically about a plane intersecting the prime meridian").

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There's a nice function I knew existed somewhere for calculating the circular mean. It's in scipy https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.circmean.html . This might be helpful but it does still need to be told the bounds to know where to wrap around so I don't know if it solves every problem but should simplify the offset case.

Id keep this function as general as possible. so I wouldn't raise an error, at worst a warning.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Have split this calculation out into its own function which now uses circmean from scipy.


# Copy cube metadata and coordinates into a new cube.
# Create a new coordinate for the coordinate along which the
Expand All @@ -252,8 +265,7 @@ def create_difference_cube(
diff_cube.add_aux_coord(coord.copy(), dims)
return diff_cube

@staticmethod
def calculate_difference(cube: Cube, coord_name: str) -> ndarray:
def calculate_difference(self, cube: Cube, coord_name: str) -> ndarray:
"""
Calculate the difference along the axis specified by the
coordinate.
Expand All @@ -268,8 +280,16 @@ def calculate_difference(cube: Cube, coord_name: str) -> ndarray:
Array after the differences have been calculated along the
specified axis.
"""
diff_axis = cube.coord_dims(coord_name)[0]
diff_along_axis = np.diff(cube.data, axis=diff_axis)
diff_axis = cube.coord(name_or_coord=coord_name)
diff_axis_number = cube.coord_dims(coord_name)[0]
diff_along_axis = np.diff(cube.data, axis=diff_axis_number)
if self._axis_wraps_around_meridian(diff_axis, cube):
first_column = cube.data[:, :1]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be clearer to have the index as 0 rather than :1.

As a separate issue it is would be better not to assume the order of the coordinates in a cube instead getting the index of the diff_axis from the cube.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Have replaced the :1/:-1 with reshapes.

I'm not sure I understand the second point. Would be good to discuss.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Have now fixed this. Does a single reshape if needed with [-1, 1] and takes the array axis to diff along from the cube dim-coord number rather than assuming the order of the coordinates on cube.data.

last_column = cube.data[:, -1:]
wrap_around_diff = np.diff(
np.hstack([last_column, first_column]), axis=diff_axis_number
)
diff_along_axis = np.hstack([diff_along_axis, wrap_around_diff])
return diff_along_axis

def process(self, cube: Cube) -> Tuple[Cube, Cube]:
Expand All @@ -291,9 +311,8 @@ def process(self, cube: Cube) -> Tuple[Cube, Cube]:
diffs = []
for axis in ["x", "y"]:
coord_name = cube.coord(axis=axis).name()
diff_cube = self.create_difference_cube(
cube, coord_name, self.calculate_difference(cube, coord_name)
)
difference = self.calculate_difference(cube, coord_name)
diff_cube = self.create_difference_cube(cube, coord_name, difference)
self._update_metadata(diff_cube, coord_name, cube.name())
diffs.append(diff_cube)
return tuple(diffs)
Expand Down
46 changes: 40 additions & 6 deletions improver_tests/utilities/test_DifferenceBetweenAdjacentGridSquares.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,24 +57,45 @@ def setUp(self):
def test_y_dimension(self):
"""Test differences calculated along the y dimension."""
points = self.cube.coord(axis="y").points
expected_y = (points[1:] + points[:-1]) / 2
expected_y_coords = (points[1:] + points[:-1]) / 2
result = self.plugin.create_difference_cube(
self.cube, "projection_y_coordinate", self.diff_in_y_array
)
self.assertIsInstance(result, Cube)
self.assertArrayAlmostEqual(result.coord(axis="y").points, expected_y)
self.assertArrayAlmostEqual(result.coord(axis="y").points, expected_y_coords)
self.assertArrayEqual(result.data, self.diff_in_y_array)

def test_x_dimension(self):
"""Test differences calculated along the x dimension."""
diff_array = np.array([[1, 1], [2, 2], [5, 5]])
points = self.cube.coord(axis="x").points
expected_x = (points[1:] + points[:-1]) / 2
expected_x_coords = (points[1:] + points[:-1]) / 2
result = self.plugin.create_difference_cube(
self.cube, "projection_x_coordinate", diff_array
)
self.assertIsInstance(result, Cube)
self.assertArrayAlmostEqual(result.coord(axis="x").points, expected_x)
self.assertArrayAlmostEqual(result.coord(axis="x").points, expected_x_coords)
self.assertArrayEqual(result.data, diff_array)

def test_x_dimension_for_circular_latlon_cube(self):
"""Test differences calculated along the x dimension."""
test_cube_data = np.array([[1, 2, 3], [2, 4, 6], [5, 10, 15]])
test_cube_x_grid_spacing = 120
test_cube = set_up_variable_cube(
test_cube_data,
"wind_speed",
"m s-1",
"latlon",
x_grid_spacing=test_cube_x_grid_spacing,
)
test_cube.coord(axis="x").circular = True
diff_array = np.array([[1, 1, -2], [2, 2, -4], [5, 5, -10]])
expected_x_coords = np.array(
[-60, 60, 180]
) # Original data are at [-120, 0, 120], therefore differences are at [-60, 60, 180].
result = self.plugin.create_difference_cube(test_cube, "longitude", diff_array)
self.assertIsInstance(result, Cube)
self.assertArrayAlmostEqual(result.coord(axis="x").points, expected_x_coords)
self.assertArrayEqual(result.data, diff_array)

def test_othercoords(self):
Expand All @@ -89,7 +110,6 @@ def test_othercoords(self):


class Test_calculate_difference(IrisTest):

"""Test the calculate_difference method."""

def setUp(self):
Expand All @@ -107,6 +127,16 @@ def test_x_dimension(self):
self.assertIsInstance(result, np.ndarray)
self.assertArrayEqual(result, expected)

def test_x_dimension_wraps_around_meridian(self):
"""Test differences calculated along the x dimension for a cube which is circular in x."""
self.cube.coord(axis="x").circular = True
expected = np.array([[1, 1, -2], [2, 2, -4], [5, 5, -10]])
result = self.plugin.calculate_difference(
self.cube, self.cube.coord(axis="x").name()
)
self.assertIsInstance(result, np.ndarray)
self.assertArrayEqual(result, expected)

def test_y_dimension(self):
"""Test differences calculated along the y dimension."""
expected = np.array([[1, 2, 3], [3, 6, 9]])
Expand Down Expand Up @@ -143,7 +173,6 @@ def test_masked_data(self):


class Test_process(IrisTest):

"""Test the process method."""

def setUp(self):
Expand Down Expand Up @@ -199,6 +228,11 @@ def test_3d_cube(self):
self.assertIsInstance(result[1], iris.cube.Cube)
self.assertArrayEqual(result[1].data, expected_y)

def test_circular_non_geographic_cube_raises_approprate_exception(self):
self.cube.coord(axis="x").circular = True
with self.assertRaises(ValueError):
self.plugin.process(self.cube)


if __name__ == "__main__":
unittest.main()
Empty file added improver_tests/utilities/test_DistanceBetweenGridPoints.py
View file
Open in desktop
Empty file.
Loading