SciTools · pelson · Jun 24, 2017 · Jun 22, 2017 · Jun 22, 2017 · Jun 22, 2017
diff --git a/lib/iris/analysis/_regrid.py b/lib/iris/analysis/_regrid.py
@@ -31,6 +31,7 @@
                                           get_xy_dim_coords, snapshot_grid)
 from iris.analysis._scipy_interpolate import _RegularGridInterpolator
 import iris.cube
+from iris.util import _meshgrid
 
 
 class RectilinearRegridder(object):
@@ -124,7 +125,7 @@ def _sample_grid(src_coord_system, grid_x_coord, grid_y_coord):
             arrays.
 
         """
-        grid_x, grid_y = np.meshgrid(grid_x_coord.points, grid_y_coord.points)
+        grid_x, grid_y = _meshgrid(grid_x_coord.points, grid_y_coord.points)
         # Skip the CRS transform if we can to avoid precision problems.
         if src_coord_system == grid_x_coord.coord_system:
             sample_grid_x = grid_x

diff --git a/lib/iris/analysis/cartography.py b/lib/iris/analysis/cartography.py
@@ -36,6 +36,7 @@
 import iris.coords
 import iris.coord_systems
 import iris.exceptions
+from iris.util import _meshgrid
 
 
 # This value is used as a fall-back if the cube does not define the earth
@@ -259,7 +260,7 @@ def get_xy_grids(cube):
 
     if x.ndim == y.ndim == 1:
         # Convert to 2D.
-        x, y = np.meshgrid(x, y)
+        x, y = _meshgrid(x, y)
     elif x.ndim == y.ndim == 2:
         # They are already in the correct shape.
         pass
@@ -284,7 +285,7 @@ def get_xy_contiguous_bounded_grids(cube):
 
     x = x_coord.contiguous_bounds()
     y = y_coord.contiguous_bounds()
-    x, y = np.meshgrid(x, y)
+    x, y = _meshgrid(x, y)
 
     return (x, y)
 
@@ -608,7 +609,7 @@ def project(cube, target_proj, nx=None, ny=None):
     source_x = lon_coord.points
     source_y = lat_coord.points
     if source_x.ndim != 2 or source_y.ndim != 2:
-        source_x, source_y = np.meshgrid(source_x, source_y)
+        source_x, source_y = _meshgrid(source_x, source_y)
 
     # Calculate target grid
     target_cs = None
@@ -1062,7 +1063,7 @@ def rotate_winds(u_cube, v_cube, target_cs):
 
     # Convert points to 2D, if not already, and determine dims.
     if x.ndim == y.ndim == 1:
-        x, y = np.meshgrid(x, y)
+        x, y = _meshgrid(x, y)
         dims = (y_dims[0], x_dims[0])
     else:
         dims = x_dims

diff --git a/lib/iris/analysis/trajectory.py b/lib/iris/analysis/trajectory.py
@@ -36,6 +36,7 @@
 from iris.analysis._interpolate_private import \
     _nearest_neighbour_indices_ndcoords, linear as linear_regrid
 from iris.analysis._interpolation import snapshot_grid
+from iris.util import _meshgrid
 
 
 class _Segment(object):
@@ -504,7 +505,7 @@ def __init__(self, src_cube, target_grid_cube):
         self.tgt_grid_shape = tgt_y_coord.shape + tgt_x_coord.shape
 
         # Calculate sample points as 2d arrays, like broadcast (NY,1)*(1,NX).
-        x_2d, y_2d = np.meshgrid(tgt_x_coord.points, tgt_y_coord.points)
+        x_2d, y_2d = _meshgrid(tgt_x_coord.points, tgt_y_coord.points)
         # Cast as a "trajectory", to suit the method used.
         self.trajectory = ((tgt_x_coord.name(), x_2d.flatten()),
                            (tgt_y_coord.name(), y_2d.flatten()))

diff --git a/lib/iris/coords.py b/lib/iris/coords.py
@@ -715,7 +715,7 @@ def _str_dates(self, dates_as_numbers):
         date_obj_array = self.units.num2date(dates_as_numbers)
         kwargs = {'separator': ', ', 'prefix': '      '}
         return np.core.arrayprint.array2string(date_obj_array,
-                                               formatter={'numpystr': str},
+                                               formatter={'all': str},
                                                **kwargs)
 
     def __str__(self):

diff --git a/lib/iris/experimental/regrid.py b/lib/iris/experimental/regrid.py
@@ -42,7 +42,7 @@
 from iris.analysis._regrid import RectilinearRegridder
 import iris.coord_systems
 import iris.cube
-from iris.util import promote_aux_coord_to_dim_coord
+from iris.util import _meshgrid, promote_aux_coord_to_dim_coord
 
 
 _Version = namedtuple('Version', ('major', 'minor', 'micro'))
@@ -755,7 +755,7 @@ def regrid_area_weighted_rectilinear_src_and_grid(src_cube, grid_cube,
 
     # Wrap up the data as a Cube.
     # Create 2d meshgrids as required by _create_cube func.
-    meshgrid_x, meshgrid_y = np.meshgrid(grid_x.points, grid_y.points)
+    meshgrid_x, meshgrid_y = _meshgrid(grid_x.points, grid_y.points)
     regrid_callback = RectilinearRegridder._regrid
     new_cube = RectilinearRegridder._create_cube(new_data, src_cube,
                                                  src_x_dim, src_y_dim,
@@ -1421,7 +1421,7 @@ def _regrid(src_data, xy_dim, src_x_coord, src_y_coord,
             src_projection, src_x_coord.points, src_y_coord.points)
 
         tgt_projection = tgt_x_coord.coord_system.as_cartopy_projection()
-        tgt_x, tgt_y = np.meshgrid(tgt_x_coord.points, tgt_y_coord.points)
+        tgt_x, tgt_y = _meshgrid(tgt_x_coord.points, tgt_y_coord.points)
         projected_tgt_grid = projection.transform_points(
             tgt_projection, tgt_x, tgt_y)
 

diff --git a/lib/iris/experimental/regrid_conservative.py b/lib/iris/experimental/regrid_conservative.py
@@ -1,4 +1,4 @@
-# (C) British Crown Copyright 2013 - 2016, Met Office
+# (C) British Crown Copyright 2013 - 2017, Met Office
 #
 # This file is part of Iris.
 #
@@ -25,9 +25,10 @@
 import cartopy.crs as ccrs
 import numpy as np
 
-from iris.analysis._interpolation import get_xy_dim_coords
 import iris
+from iris.analysis._interpolation import get_xy_dim_coords
 from iris.analysis._regrid import RectilinearRegridder
+from iris.util import _meshgrid
 
 
 #: A static Cartopy Geodetic() instance for transforming to true-lat-lons.
@@ -81,8 +82,8 @@ def _make_esmpy_field(x_coord, y_coord, ref_name='field',
     grid = ESMF.Grid(dims)
 
     # Get all cell corner coordinates as true-lat-lons
-    x_bounds, y_bounds = np.meshgrid(x_coord.contiguous_bounds(),
-                                     y_coord.contiguous_bounds())
+    x_bounds, y_bounds = _meshgrid(x_coord.contiguous_bounds(),
+                                   y_coord.contiguous_bounds())
     grid_crs = x_coord.coord_system.as_cartopy_crs()
     lon_bounds, lat_bounds = _convert_latlons(grid_crs, x_bounds, y_bounds)
 
@@ -113,7 +114,7 @@ def _make_esmpy_field(x_coord, y_coord, ref_name='field',
     x_centres = 0.5 * (x_centres[:-1] + x_centres[1:])
     y_centres = y_coord.contiguous_bounds()
     y_centres = 0.5 * (y_centres[:-1] + y_centres[1:])
-    x_points, y_points = np.meshgrid(x_centres, y_centres)
+    x_points, y_points = _meshgrid(x_centres, y_centres)
     lon_points, lat_points = _convert_latlons(grid_crs, x_points, y_points)
 
     # Add grid 'coord' element for centres + fill with centre-points values.

diff --git a/lib/iris/plot.py b/lib/iris/plot.py
@@ -49,6 +49,7 @@
 from iris.exceptions import IrisError
 # Importing iris.palette to register the brewer palettes.
 import iris.palette
+from iris.util import _meshgrid
 
 
 # Cynthia Brewer citation text.
@@ -703,16 +704,16 @@ def _map_common(draw_method_name, arg_func, mode, cube, plot_defn,
     y_coord, x_coord = plot_defn.coords
     if mode == iris.coords.POINT_MODE:
         if x_coord.ndim == y_coord.ndim == 1:
-            x, y = np.meshgrid(x_coord.points, y_coord.points)
+            x, y = _meshgrid(x_coord.points, y_coord.points)
         elif x_coord.ndim == y_coord.ndim == 2:
             x = x_coord.points
             y = y_coord.points
         else:
             raise ValueError("Expected 1D or 2D XY coords")
     else:
         try:
-            x, y = np.meshgrid(x_coord.contiguous_bounds(),
-                               y_coord.contiguous_bounds())
+            x, y = _meshgrid(x_coord.contiguous_bounds(),
+                             y_coord.contiguous_bounds())
         # Exception translation.
         except iris.exceptions.CoordinateMultiDimError:
             raise ValueError("Could not get XY grid from bounds. "

diff --git a/lib/iris/tests/results/analysis/rotated_pole.0.rlat.json b/lib/iris/tests/results/analysis/rotated_pole.0.rlat.json
@@ -1 +1 @@
-{"std": 3.7195861803241614, "min": -5.593200206756592, "max": 8.72130012512207, "shape": [93, 75], "masked": false, "mean": 1.6348200228305594}
+{"std": 3.7195863723754883, "min": -5.5932002067565918, "max": 8.7213001251220703, "shape": [93, 75], "masked": false, "mean": 1.6348202228546143}
diff --git a/lib/iris/tests/unit/data_manager/test_DataManager.py b/lib/iris/tests/unit/data_manager/test_DataManager.py
@@ -662,7 +662,7 @@ def test_with_lazy_mask_array__masked(self):
         self.assertIsInstance(result, ma.MaskedArray)
         self.assertIsNone(dm._realised_dtype)
         self.assertEqual(dm.dtype, self.realised_dtype)
-        self.assertArrayEqual(result, self.lazy_mask_array_masked.compute())
+        self.assertArrayEqual(result, self.mask_array_masked)
 
     def test_with_real_masked_constant(self):
         masked_data = ma.masked_array([666], mask=True, dtype=np.dtype('f8'))

diff --git a/lib/iris/util.py b/lib/iris/util.py
@@ -1657,3 +1657,23 @@ def demote_dim_coord_to_aux_coord(cube, name_or_coord):
     cube.remove_coord(dim_coord)
 
     cube.add_aux_coord(dim_coord, coord_dim)
+
+
+@functools.wraps(np.meshgrid)
+def _meshgrid(*xi, **kwargs):
+    """
+    @numpy v1.13, the dtype of each output nD coordinate is the same as its
+    associated input 1D coordinate. This is not the case prior to numpy v1.13,
+    where the output dtype is cast up to its highest resolution, regardlessly.
+
+    This convenience function ensures consistent meshgrid behaviour across
+    numpy versions.
+
+    Reference: https://github.com/numpy/numpy/pull/5302
+
+    """
+    mxi = np.meshgrid(*xi, **kwargs)
+    for i, (mxii, xii) in enumerate(zip(mxi, xi)):
+        if mxii.dtype != xii.dtype:
+            mxi[i] = mxii.astype(xii.dtype)
+    return mxi
Original file line number	Diff line number	Diff line change
		@@ -1 +1 @@
		{"std": 3.7195861803241614, "min": -5.593200206756592, "max": 8.72130012512207, "shape": [93, 75], "masked": false, "mean": 1.6348200228305594}
		{"std": 3.7195863723754883, "min": -5.5932002067565918, "max": 8.7213001251220703, "shape": [93, 75], "masked": false, "mean": 1.6348202228546143}
Copy link Member pp-mo Jun 23, 2017 • edited Loading Choose a reason for hiding this comment The reason will be displayed to describe this comment to others. Learn more. It might have been "nicer" here to reduce the tolerance on the test for when it happens again (!). I see these values only just fail the test at a relative tolerance of 1e-7, as default for `assertArrayAllClose` (+hence `assertDataAlmostEqual`) : mean value relative change is approx 1.22e-7. You could maybe add something like `, reltol=5.0e-6` to the 4 calls here You could do that instead, or as well if you prefer (changes are very small after all). Just a suggestion -- skip if you don't want. Copy link Member Author bjlittle Jun 23, 2017 Choose a reason for hiding this comment The reason will be displayed to describe this comment to others. Learn more. @pp-mo I'm really against fudging the tolerances. We took that step previously with testing and it ultimately lead to a bad place, particularly with graphical testing. So, I'd rather not take that option, unless you consider it unwise. Copy link Member pp-mo Jun 23, 2017 Choose a reason for hiding this comment The reason will be displayed to describe this comment to others. Learn more. I don't consider it fudging as long as the value used remains appropriate. I know that was not always the case for our image comparison tolerance, where it got stupidly huge for a while... ! The default value of 1e-7 is designed to cope with float32 rounding errors, but can easily be a bit fine for testing other calculated values, like coordinate projections. For instance, I just encountered some other tests failing at that level when I was messing about with different dask chunkings.