Make difference module handle circular cubes #1990
Conversation
…ted correctly but need to add extra dimension for the wrap-around distance.
MO-PeterJordan
changed the title
improver/utilities/spatial.py
Outdated
| extra_mean_point = ( | ||
| np.mean([points[-1], (points[0] + max_value)]) % max_value | ||
| ) | ||
| mean_points = np.hstack([mean_points, extra_mean_point]) |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
...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").
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
Have split this calculation out into its own function which now uses circmean from scipy.
MO-PeterJordan
changed the title
MO-PeterJordan
changed the title
…actions despite running locally and pushing
improver/utilities/spatial.py
Outdated
| extra_mean_point = ( | ||
| np.mean([points[-1], (points[0] + max_value)]) % max_value | ||
| ) | ||
| mean_points = np.hstack([mean_points, extra_mean_point]) |
There was a problem hiding this comment.
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.
improver/utilities/spatial.py
Outdated
| 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] |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
Have replaced the :1/:-1 with reshapes.
I'm not sure I understand the second point. Would be good to discuss.
There was a problem hiding this comment.
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.
improver_tests/utilities/test_DifferenceBetweenAdjacentGridSquares.py
Outdated
Show resolved
Hide resolved
improver/utilities/spatial.py
Outdated
| extra_mean_point = ( | ||
| np.mean([points[-1], (points[0] + max_value)]) % max_value | ||
| ) | ||
| mean_points = np.hstack([mean_points, extra_mean_point]) |
There was a problem hiding this comment.
Have split this calculation out into its own function which now uses circmean from scipy.
improver/utilities/spatial.py
Outdated
| 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] |
There was a problem hiding this comment.
Have replaced the :1/:-1 with reshapes.
I'm not sure I understand the second point. Would be good to discuss.
improver_tests/utilities/test_DifferenceBetweenAdjacentGridSquares.py
Outdated
Show resolved
Hide resolved
improver/utilities/spatial.py
Outdated
| first_column = cube.data[:, 0].reshape([-1, 1]) | ||
| last_column = cube.data[:, -1].reshape([-1, 1]) | ||
| wrap_around_diff = np.diff( |
There was a problem hiding this comment.
This assumes that the cube is set up with x as rows and y as columns but this may not be true.
Can use an Improver module to force the cube to be x as rows and y as columns but have to make sure I've put it back again. Alternatively Might be able to work out correct index here on the fly.
There was a problem hiding this comment.
Have added a test for a cube using the opposite lat/long axis choices and fixed this.
Below is just my chain of thought/sanity check. Feel free to ignore.
...Wait... I produced the flipped code by transposing the cube. This flips which axis of cube.data each dim coord is assigned to, and also transposes the numpy array, so that the data represented by the cube is the same, but the longitudes now increase along the columns and the latitudes along the rows... which is exactly what I want. Okay, good.
…umns not rows and vice versa for latitude). Fixed code to handle this case.
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #1990 +/- ##
==========================================
+ Coverage 98.39% 98.41% +0.01%
==========================================
Files 124 128 +4
Lines 12212 12482 +270
==========================================
+ Hits 12016 12284 +268
- Misses 196 198 +2 ☔ View full report in Codecov by Sentry. |
|
Superseded. |
Split out from #1985 to support calculating gradients for cubes which wrap around the meridian.
Testing:
CLA