Gridcell rotations #3096

pp-mo · 2018-07-26T13:09:45Z

**WIP - DO NOT MERGE probably can now merge, once caveats are established
Here is a first look at the gridcell angles code.

The tests are currently a tangled mess of exploratory code + not fit for human consumption !
Eventually, I will have just the 2 new routines and proper unit tests for them.

There is a fair bit that still needs to be done, even in the main code :

improve documentation generally
provide for coords in alternative coord-systems, by projecting to ordinary lat-lon
review how coords are identified (axis rather than 'longitude' + 'latitude' names)

…phy, with tests: INCOMPLETE WIP

stickler-ci · 2018-07-26T13:10:12Z

lib/iris/analysis/_grid_angles.py

+    if old_style:
+        mid_lons = np.deg2rad(q[0])
+
+        pr =  _3d_xyz_from_latlon(r[0], r[1]) - _3d_xyz_from_latlon(p[0], p[1])


E222 multiple spaces after operator

stickler-ci · 2018-07-26T13:10:13Z

lib/iris/analysis/_grid_angles.py

+                                                 lmb_hatvec_y,
+                                                 lmb_hatvec_z)])
+
+        pr =  _3d_xyz_from_latlon(r[0], r[1]) - _3d_xyz_from_latlon(p[0], p[1])


E222 multiple spaces after operator

stickler-ci · 2018-07-26T13:10:13Z

lib/iris/analysis/cartography.py

@@ -37,6 +37,9 @@
 import iris.coord_systems
 import iris.exceptions
 from iris.util import _meshgrid
+from ._grid_angles import (


F401 '._grid_angles.true_vectors_from_grid_vectors as rotate_grid_vectors' imported but unused
F401 '._grid_angles.gridcell_angles' imported but unused

stickler-ci · 2018-07-26T13:10:13Z

lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py

+import iris.tests as tests
+
+import numpy as np
+import numpy.ma as ma


F401 'numpy.ma' imported but unused

stickler-ci · 2018-07-26T13:10:13Z

lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py

+
+import cartopy.crs as ccrs
+from iris.cube import Cube
+from iris.coords import DimCoord, AuxCoord


F401 'iris.coords.DimCoord' imported but unused

stickler-ci · 2018-07-26T13:10:16Z

lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py

+##        for iy in range(ny):
+##            for ix in range(nx):
+##                plt.plot([x0s[iy, ix], x1s[iy, ix]],
+##                         [y0s[iy, ix], y1s[iy, ix]],


E266 too many leading '#' for block comment

stickler-ci · 2018-07-26T13:10:16Z

lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py

+##            for ix in range(nx):
+##                plt.plot([x0s[iy, ix], x1s[iy, ix]],
+##                         [y0s[iy, ix], y1s[iy, ix]],
+##                         'o-', markersize=4., markeredgewidth=0.,


E266 too many leading '#' for block comment

stickler-ci · 2018-07-26T13:10:16Z

lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py

+##                plt.plot([x0s[iy, ix], x1s[iy, ix]],
+##                         [y0s[iy, ix], y1s[iy, ix]],
+##                         'o-', markersize=4., markeredgewidth=0.,
+##                         color='green', # scale_units='xy', scale=scale,


E266 too many leading '#' for block comment

stickler-ci · 2018-07-26T13:10:16Z

lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py

+##                         [y0s[iy, ix], y1s[iy, ix]],
+##                         'o-', markersize=4., markeredgewidth=0.,
+##                         color='green', # scale_units='xy', scale=scale,
+##                         transform=data_proj)


E266 too many leading '#' for block comment

stickler-ci · 2018-07-26T13:10:16Z

lib/iris/tests/unit/analysis/cartography/test_gridcell_angles.py

+
+
+
+        ax.set_global()


E303 too many blank lines (3)

lbdreyer

Not done! I'm making my way though the _angle function and then just have tests to look through, but I thought I'd submit this for now

lbdreyer · 2018-07-27T12:06:23Z

lib/iris/analysis/_grid_angles.py

+    uu, vv = mags * np.cos(angs), mags * np.sin(angs)
+
+    # Promote all to masked arrays, and also apply mask at bad (NaN) angles.
+    mask = np.isnan(aa)


Can aa i.e. the grid cell angles also be masked, or just nans?

I hadn't thought of that.
I think it doesn't much matter, as the results would just be NaN at those points.

... however, now I'm wondering if we should return masked arrays in every case.
Iris is usually more intelligent than that.

lbdreyer · 2018-07-27T12:24:19Z

lib/iris/analysis/_grid_angles.py

+        Valid values are 'lower-left, lower-right', which takes the angle from
+        the lower left to the lower right corner, and 'mid-lhs, mid-rhs' which
+        takes an angles between the average of the left-hand and right-hand
+        pairs of corners.  The default is 'mid-lhs, mid-rhs'.


Any chance it would ever need to be upper-left, upper right?

The original code I was working from (ex-IDL) used only point values.
I extended it to use bounds because that is much more logical, and also works properly for regional data.
What I found was that, while 'mid-lhs..mid-rhs' was a more obvious choice, 'bottom-left..bottom-right' was closer to the 'official' cell angles in the model output file (when included).
So I hedged my bets by making it switchable, but only supported 2 options = 'seems most logical' and 'best match to official values'.

We can perhaps standardise on just the 'mid-lhs..mid-rhs' option, if I could convince myself that the differences are not significant. I don't really want to make 'bottom-left..bottom-right' the default, as it seems the poorer choice.

standardise on just the 'mid-lhs..mid-rhs' option, if I could convince myself that the differences are not significant

But unfortunately I checked, and the discrepancies are significant in some areas. At least a few degrees for practical purposes.
Also, the errors for quasi-ideal square-ish gridcells rotated at various angles seem to be ~linear in the size of the cells for the 'bottom-left..bottom-right' method while ~square for the 'mid-lhs..mid-rhs'.
So that means I'm not comfortable 'standardising' on either of these exclusively 😞

lbdreyer · 2018-07-27T12:29:45Z

lib/iris/analysis/_grid_angles.py

+
+    """
+    cube = None
+    if hasattr(x, 'add_aux_coord'):


Why aren't you going isinstance?

Because duck typing is "Better" ™️ !

lbdreyer · 2018-07-27T12:30:46Z

lib/iris/analysis/_grid_angles.py

+        raise ValueError(msg.format(type(x), type(y)))
+
+    x_coord, y_coord = None, None
+    if hasattr(x, 'bounds') and hasattr(y, 'bounds'):


As before, why not check for isinstance to check it's a cube?

lbdreyer · 2018-07-27T12:40:29Z

lib/iris/analysis/_grid_angles.py

+        msg = 'Input {!r} is a Coordinate, but {!r} is not.'
+        raise ValueError(*is_and_not)
+
+    # Now have either 2 points arrays or 2 bounds arrays.


Can this be more explicit, i.e.
# Now we should have either 2 points arrays of shape (nx, ny) or 2 bounds arrays of shape (nx, ny, 4)

lbdreyer · 2018-07-27T12:45:00Z

lib/iris/analysis/_grid_angles.py

+    if x.ndim == 2:
+        # PROBLEM: we can't use this if data is not full-longitudes,
+        # i.e. rhs of array must connect to lhs (aka 'circular' coordinate).
+        # But we have no means of checking that ?


⚠️
This is concerning

Yes, but that's how it is.
I don't see any practical alternative.
Instead, I have stated this behaviour in the docstring.

The only extra thing I can imagine to do here would be to insist on 'circular' if we are passed an x-coord, but I'm very down on that : As a rare example of non-CF metadata, I think 'circular' it is a bit of a mistake, and I don't think Iris uses the concept at all consistently -- I guess it is probably is not set in many cases where it would be appropriate.

lbdreyer · 2018-07-27T13:10:53Z

lib/iris/analysis/_grid_angles.py

+        # NOTE: so with bounds, we *don't* need full circular longitudes.
+        xyz = _3d_xyz_from_latlon(x, y)
+        angle_boundpoints_vals = {'mid-lhs, mid-rhs': '03_to_12',
+                                  'lower-left, lower-right': '0_to_1'}


Kinda confusing names 03_to_12 and 0_to_1. I initially thought it may have something to do with a clock.
I understand them now, but I will try to think up something informative

lbdreyer · 2018-07-27T13:19:53Z

lib/iris/analysis/_grid_angles.py

+                   'must be one of {}')
+            raise ValueError(msg.format(cell_angle_boundpoints,
+                                        list(angle_boundpoints_vals.keys())))
+        if not x_coord:


Maybe better to do
if not x_coord and y_coord
I'm not sure if you catch the example someone supplying giving gridcell_angles a coord and an array as arguments. So maybe it would, instead, be better to check that earlier?

The code above checked for "one coord + one non-coord", and will error that, here.

lbdreyer · 2018-07-27T13:28:10Z

lib/iris/analysis/_grid_angles.py

+
+        angles : (2-dimensional cube)
+
+            Cube of angles of grid-x vector from true Eastward direction for


How about, if the input was a cube, the grid cell angles are an aux coord on the original cube? Or is that too confusing?

It just seems wasteful to have a copy of the lat and lon coord somewhere else. I guess this is an example of where shared coordinates would be useful

We can have gridcell angles in the file, but as they are not a CF concept they never appear as an aux-coord, but instead just as a separate cube.

Talking of inefficient, of course the "angles" cube has another copy of the lats + lons with all their bounds!
We really do need to fix this ...

lbdreyer · 2018-07-27T13:48:09Z

lib/iris/analysis/_grid_angles.py

+import iris
+
+
+def _3d_xyz_from_latlon(lon, lat):


These functions ( _3d_xyz_from_latlon and '_latlon_from_xyz') seems quite general purpose functions (not just specific to the case of rotating grid cells)... could we not put them somewhere else so they can be used by other things?

@pp-mo have you logged this suggestion in a labelled issue yet? Might be better to merge first and think about where things should live later.

pp-mo · 2018-07-30T10:26:33Z

**WIP - DO NOT MERGE probably can now merge, once caveats are established

TODO:

I have already prepared somemore changes + tidying to this
- these are pending in another branch : Vector plots 2 pp-mo/iris#45
- note: just because ... I don't want to mess up this discussion yet
once we have resolved the outstanding stuff here, I will merge that in
also various additional (minor) testing is wanted
which I propose to defer for now

lbdreyer · 2018-07-30T12:22:38Z

lib/iris/analysis/_grid_angles.py

+    Return angle (in _radians_) of grid wrt local east.
+    Anticlockwise +ve, as usual.
+    {P, Q, R} are consecutive points in the same row,
+    eg {v(i,j),f(i,j),v(i+1,j)}, or {T(i-1,j),T(i,j),T(i+1,j)}


Should this not be:
{v(i,j),v(i,j),v(i+1,j)}

Why are v and T used here?

Also, I assume these function docs strings are not quite finalised yet? So I will reserve comments on the wording till then

This is inherited docs + it does need improving.
Let's defer that : #3109

( I think the 'v', 'f', and 'T' are referring to different types of model value sited at different gridcell locations (i.e. stagger), so 'v' is on a v-grid, while 'T' (and 'f'?) will be on the 'T' grid. )

lbdreyer · 2018-07-30T12:35:04Z

lib/iris/analysis/_grid_angles.py

+    """
+#    old_style = True
+    old_style = False
+    if old_style:


How do you plan to determine whether or not old_style should be used?

Sorry, this was an alternative calculation that I used for evaluation.
I have since removed it (cf. pp-mo#45)

Though... I accept that I still don't have a full understanding of this calculation 😬 .
We should probably make more effort to explain exactly what this is doing.
Put that also down to #3109 !

corinnebosley · 2018-07-31T11:02:10Z

lib/iris/analysis/_grid_angles.py

+    Args:
+
+    * xyz: (array)
+        positions array, of dims (3, <others>), where index 0 maps x/y/z.


This is a bit confusing. I know that the docs in general are not finished though, so I won't go on.

corinnebosley · 2018-07-31T11:02:27Z

lib/iris/analysis/_grid_angles.py

+    Returns:
+
+        lonlat : (array)
+            spherical angles, of dims (2, <others>), in radians.


so many commas...

We need to work out how to capture outstanding review issues on things we merged in advance.
Perhaps like this ... #3114 ?

pp-mo added 2 commits July 23, 2018 16:10

Add gridcell_angles and rotate_grid_vectors to iris.analysis.cartogra…

52496b9

…phy, with tests: INCOMPLETE WIP

Roughly working, snapshotted with complex test plot code, to be reduced.

2e0da2e

stickler-ci reviewed Jul 26, 2018

View reviewed changes

Small improvements.

6b9bfbb

pp-mo force-pushed the vector_plots branch from 53e0892 to 6b9bfbb Compare July 26, 2018 15:16

lbdreyer reviewed Jul 27, 2018

View reviewed changes

pp-mo changed the base branch from master to 2d_coords July 30, 2018 10:10

pp-mo added the 2d_coords label Jul 30, 2018

pp-mo mentioned this pull request Jul 30, 2018

Vector plots 2 pp-mo/iris#45

Closed

lbdreyer reviewed Jul 30, 2018

View reviewed changes

pp-mo mentioned this pull request Jul 30, 2018

Improve comments describing calculation in iris.analysis.cartography.grid_angles #3109

Closed

corinnebosley reviewed Jul 31, 2018

View reviewed changes

pp-mo mentioned this pull request Jul 31, 2018

Outstanding review issues from #3096 (gridcell rotations) #3114

Closed

lbdreyer merged commit c098c5a into SciTools:2d_coords Jul 31, 2018

pp-mo mentioned this pull request Jul 31, 2018

Vector plots 2 #3120

Merged

pp-mo deleted the vector_plots branch January 8, 2019 12:57


		angles : (2-dimensional cube)

		Cube of angles of grid-x vector from true Eastward direction for

Gridcell rotations #3096

Gridcell rotations #3096

Uh oh!

Conversation

pp-mo commented Jul 26, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

lbdreyer left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pp-mo Jul 30, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pp-mo Jul 27, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pp-mo Jul 27, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pp-mo commented Jul 30, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

pp-mo commented Jul 26, 2018 •

edited

Loading

pp-mo Jul 30, 2018 •

edited

Loading

pp-mo Jul 27, 2018 •

edited

Loading

pp-mo Jul 27, 2018 •

edited

Loading

pp-mo commented Jul 30, 2018 •

edited

Loading

pp-mo Jul 30, 2018 •

edited

Loading