Port the changes from iris.fileformats.grib back into iris_grib.

iris_grib/__init__.py

@@ -29,14 +29,14 @@
 import math  # for fmod
 import warnings
 
-import biggus
 import cartopy.crs as ccrs
 import cf_units
 import gribapi
 import numpy as np
 import numpy.ma as ma
 
 import iris
+from iris._lazy_data import as_lazy_data, convert_nans_array
 import iris.coord_systems as coord_systems
 from iris.exceptions import TranslationError, NotYetImplementedError
 
@@ -99,12 +99,11 @@
 class GribDataProxy(object):
     """A reference to the data payload of a single Grib message."""
 
-    __slots__ = ('shape', 'dtype', 'fill_value', 'path', 'offset')
+    __slots__ = ('shape', 'dtype', 'path', 'offset')
 
-    def __init__(self, shape, dtype, fill_value, path, offset):
+    def __init__(self, shape, dtype, path, offset):
         self.shape = shape
         self.dtype = dtype
-        self.fill_value = fill_value
         self.path = path
         self.offset = offset
 
@@ -123,7 +122,7 @@ def __getitem__(self, keys):
 
     def __repr__(self):
         msg = '<{self.__class__.__name__} shape={self.shape} ' \
-            'dtype={self.dtype!r} fill_value={self.fill_value!r} ' \
+            'dtype={self.dtype!r} ' \
             'path={self.path!r} offset={self.offset}>'
         return msg.format(self=self)
 
@@ -146,6 +145,7 @@ class GribWrapper(object):
     def __init__(self, grib_message, grib_fh=None):
         """Store the grib message and compute our extra keys."""
         self.grib_message = grib_message
+        self.realised_dtype = np.array([0.]).dtype
 
         if self.edition != 1:
             emsg = 'GRIB edition {} is not supported by {!r}.'
@@ -184,12 +184,14 @@ def __init__(self, grib_message, grib_fh=None):
             # The byte offset requires to be reset back to the first byte
             # of this message. The file pointer offset is always at the end
             # of the current message due to the grib-api reading the message.
-            proxy = GribDataProxy(shape, np.zeros(0).dtype, np.nan,
-                                  grib_fh.name,
+            proxy = GribDataProxy(shape, self.realised_dtype, grib_fh.name,
                                   offset - message_length)
-            self._data = biggus.NumpyArrayAdapter(proxy)
+            self._data = as_lazy_data(proxy)
         else:
-            self.data = _message_values(grib_message, shape)
+            values_array = _message_values(grib_message, shape)
+            # mask where the values are nan
+            self.data = convert_nans_array(values_array,
+                                           nans_replacement=ma.masked)
 
     def _confirm_in_scope(self):
         """Ensure we have a grib flavour that we choose to support."""
@@ -640,6 +642,19 @@ def _get_verification_date(self):
         # Return validity_time = (reference_time + start_offset*time_unit).
         return reference_date_time + interval_delta
 
+    @property
+    def bmdi(self):
+        # Not sure of any cases where GRIB provides a fill value.
+        # Default for fill value is None.
+        return None
+
+    def core_data(self):
+        try:
+            data = self._data
+        except AttributeError:
+            data = self.data
+        return data
+
     def phenomenon_points(self, time_unit):
         """
         Return the phenomenon time point offset from the epoch time reference
@@ -683,10 +698,6 @@ def _message_values(grib_message, shape):
     data = gribapi.grib_get_double_array(grib_message, 'values')
     data = data.reshape(shape)
 
-    # Handle missing values in a sensible way.
-    mask = np.isnan(data)
-    if mask.any():
-        data = ma.array(data, mask=mask, fill_value=np.nan)
     return data
 
 

iris_grib/_load_convert.py

@@ -414,7 +414,7 @@ def ellipsoid(shapeOfTheEarth, major, minor, radius):
 
     """
     # Supported shapeOfTheEarth values.
-    if shapeOfTheEarth not in (0, 1, 3, 6, 7):
+    if shapeOfTheEarth not in (0, 1, 2, 3, 4, 5, 6, 7):
         msg = 'Grid definition section 3 contains an unsupported ' \
             'shape of the earth [{}]'.format(shapeOfTheEarth)
         raise TranslationError(msg)
@@ -425,25 +425,36 @@ def ellipsoid(shapeOfTheEarth, major, minor, radius):
     elif shapeOfTheEarth == 1:
         # Earth assumed spherical with radius specified (in m) by
         # data producer.
-        if radius is ma.masked:
+        if ma.is_masked(radius):
             msg = 'Ellipsoid for shape of the earth {} requires a' \
                 'radius to be specified.'.format(shapeOfTheEarth)
             raise ValueError(msg)
         result = icoord_systems.GeogCS(radius)
+    elif shapeOfTheEarth == 2:
+        # Earth assumed oblate spheroid with size as determined by IAU in 1965.
+        result = icoord_systems.GeogCS(6378160, inverse_flattening=297.0)
     elif shapeOfTheEarth in [3, 7]:
         # Earth assumed oblate spheroid with major and minor axes
         # specified (in km)/(in m) by data producer.
         emsg_oblate = 'Ellipsoid for shape of the earth [{}] requires a' \
             'semi-{} axis to be specified.'
-        if major is ma.masked:
+        if ma.is_masked(major):
             raise ValueError(emsg_oblate.format(shapeOfTheEarth, 'major'))
-        if minor is ma.masked:
+        if ma.is_masked(minor):
             raise ValueError(emsg_oblate.format(shapeOfTheEarth, 'minor'))
         # Check whether to convert from km to m.
         if shapeOfTheEarth == 3:
             major *= 1000
             minor *= 1000
         result = icoord_systems.GeogCS(major, minor)
+    elif shapeOfTheEarth == 4:
+        # Earth assumed oblate spheroid as defined in IAG-GRS80 model.
+        result = icoord_systems.GeogCS(6378137,
+                                       inverse_flattening=298.257222101)
+    elif shapeOfTheEarth == 5:
+        # Earth assumed represented by WGS84 (as used by ICAO since 1998).
+        result = icoord_systems.GeogCS(6378137,
+                                       inverse_flattening=298.257223563)
     elif shapeOfTheEarth == 6:
         # Earth assumed spherical with radius of 6 371 229.0m
         result = icoord_systems.GeogCS(6371229)

iris_grib/_save_rules.py

@@ -1156,22 +1156,18 @@ def product_definition_section(cube, grib):
 
 def data_section(cube, grib):
     # Masked data?
-    if isinstance(cube.data, ma.core.MaskedArray):
-        # What missing value shall we use?
-        if not np.isnan(cube.data.fill_value):
-            # Use the data's fill value.
-            fill_value = float(cube.data.fill_value)
-        else:
-            # We can't use the data's fill value if it's NaN,
+    if ma.isMaskedArray(cube.data):
+        fill_value = cube.fill_value
+        if fill_value is None or np.isnan(cube.fill_value):
+            # We can't use the cube's fill value if it's NaN,
             # the GRIB API doesn't like it.
             # Calculate an MDI outside the data range.
             min, max = cube.data.min(), cube.data.max()
             fill_value = min - (max - min) * 0.1
-        # Prepare the unmaksed data array, using fill_value as the MDI.
-        data = cube.data.filled(fill_value)
     else:
         fill_value = None
-        data = cube.data
+
+    data = cube.data
 
     # units scaling
     grib2_info = gptx.cf_phenom_to_grib2_info(cube.standard_name,

iris_grib/message.py

@@ -1,4 +1,4 @@
-# (C) British Crown Copyright 2014 - 2016, Met Office
+# (C) British Crown Copyright 2014 - 2017, Met Office
 #
 # This file is part of iris-grib.
 #
@@ -26,10 +26,11 @@
 from collections import namedtuple
 import re
 
-import biggus
 import gribapi
 import numpy as np
+import numpy.ma as ma
 
+from iris._lazy_data import array_masked_to_nans, as_lazy_data
 from iris.exceptions import TranslationError
 
 
@@ -88,7 +89,7 @@ def __init__(self, raw_message, recreate_raw, file_ref=None):
         """
         # A RawGribMessage giving gribapi access to the original grib message.
         self._raw_message = raw_message
-        # A _MessageLocation which biggus uses to read the message data array,
+        # A _MessageLocation which dask uses to read the message data array,
         # by which time this message may be dead and the original grib file
         # closed.
         self._recreate_raw = recreate_raw
@@ -113,10 +114,23 @@ def sections(self):
         """
         return self._raw_message.sections
 
+    @property
+    def bmdi(self):
+        # Not sure of any cases where GRIB provides a fill value.
+        # Default for fill value is None.
+        return None
+
+    @property
+    def realised_dtype(self):
+        return np.dtype('f8')
+
+    def core_data(self):
+        return self.data
+
     @property
     def data(self):
         """
-        The data array from the GRIB message as a biggus Array.
+        The data array from the GRIB message as a dask Array.
 
         The shape of the array will match the logical shape of the
         message's grid. For example, a simple global grid would be
@@ -150,9 +164,9 @@ def data(self):
                 shape = (grid_section['numberOfDataPoints'],)
             else:
                 shape = (grid_section['Nj'], grid_section['Ni'])
-            proxy = _DataProxy(shape, np.dtype('f8'), np.nan,
+            proxy = _DataProxy(shape, self.realised_dtype, np.nan,
                                self._recreate_raw)
-            data = biggus.NumpyArrayAdapter(proxy)
+            data = as_lazy_data(proxy)
         else:
             fmt = 'Grid definition template {} is not supported'
             raise TranslationError(fmt.format(template))
@@ -180,12 +194,13 @@ def __call__(self):
 class _DataProxy(object):
     """A reference to the data payload of a single GRIB message."""
 
-    __slots__ = ('shape', 'dtype', 'fill_value', 'recreate_raw')
+    __slots__ = ('shape', 'dtype', 'recreate_raw')
 
     def __init__(self, shape, dtype, fill_value, recreate_raw):
+        # TODO: I (@pelson) have no idea why fill_value remains an argument as
+        # it isn't used. It was copied verbatim from iris' dask branch.
         self.shape = shape
         self.dtype = dtype
-        self.fill_value = fill_value
         self.recreate_raw = recreate_raw
 
     @property
@@ -245,20 +260,21 @@ def __getitem__(self, keys):
                 # Only the non-masked values are included in codedValues.
                 _data = np.empty(shape=bitmap.shape)
                 _data[bitmap.astype(bool)] = data
-                # `np.ma.masked_array` masks where input = 1, the opposite of
-                # the behaviour specified by the GRIB spec.
-                data = np.ma.masked_array(_data, mask=np.logical_not(bitmap))
+                # Use nan where input = 1, the opposite of the behaviour
+                # specified by the GRIB spec.
+                _data[np.logical_not(bitmap.astype(bool))] = np.nan
+                data = _data
             else:
                 msg = 'Shapes of data and bitmap do not match.'
                 raise TranslationError(msg)
 
         data = data.reshape(self.shape)
+
         return data.__getitem__(keys)
 
     def __repr__(self):
         msg = '<{self.__class__.__name__} shape={self.shape} ' \
-            'dtype={self.dtype!r} fill_value={self.fill_value!r} ' \
-            'recreate_raw={self.recreate_raw!r} '
+            'dtype={self.dtype!r} recreate_raw={self.recreate_raw!r} '
         return msg.format(self=self)
 
     def __getstate__(self):

iris_grib/tests/results/unit/load_cubes/load_cubes/reduced_raw.cml

@@ -1,6 +1,6 @@
 <?xml version="1.0" ?>
 <cubes xmlns="urn:x-iris:cubeml-0.2">
-  <cube standard_name="geopotential" units="m2 s-2">
+  <cube core-dtype="float64" dtype="float64" standard_name="geopotential" units="m2 s-2">
     <attributes>
       <attribute name="centre" value="European Centre for Medium Range Weather Forecasts"/>
     </attributes>

iris_grib/tests/unit/grib1_load_rules/test_grib1_convert.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-grib.
 #
@@ -61,7 +61,9 @@ def assert_bounded_message(self, **kwargs):
                       '_forecastTimeUnit': 'hours',
                       'phenomenon_bounds': lambda u: (80, 120),
                       '_phenomenonDateTime': -1,
-                      'table2Version': 9999}
+                      'table2Version': 9999,
+                      '_originatingCentre': 'xxx',
+                      }
         attributes.update(kwargs)
         message = mock.Mock(**attributes)
         self._test_for_coord(message, grib1_convert, self.is_forecast_period,

iris_grib/tests/unit/load_convert/test_ellipsoid.py

@@ -1,4 +1,4 @@
-# (C) British Crown Copyright 2014 - 2016, Met Office
+# (C) British Crown Copyright 2014 - 2017, Met Office
 #
 # This file is part of iris-grib.
 #
@@ -42,7 +42,7 @@
 
 class Test(tests.IrisGribTest):
     def test_shape_unsupported(self):
-        unsupported = [2, 4, 5, 8, 9, 10, MDI]
+        unsupported = [8, 9, 10, MDI]
         emsg = 'unsupported shape of the earth'
         for shape in unsupported:
             with self.assertRaisesRegexp(TranslationError, emsg):