Fix erroneous use of grid_latitude and grid_longitude and cleaned ocean grid fixes

esmvalcore/cmor/_fixes/cmip5/bcc_csm1_1.py

@@ -1,68 +1,7 @@
 """Fixes for bcc-csm1-1."""
-import numpy as np
-from scipy.interpolate import InterpolatedUnivariateSpline
-from scipy.ndimage import map_coordinates
-
-from ..common import ClFixHybridPressureCoord
-from ..fix import Fix
-
+from ..common import ClFixHybridPressureCoord, OceanFixGrid
 
 Cl = ClFixHybridPressureCoord
 
 
-class Tos(Fix):
-    """Fixes for tos."""
-
-    def fix_data(self, cube):
-        """Fix data.
-
-        Calculate missing latitude/longitude boundaries using interpolation.
-
-        Parameters
-        ----------
-        cube: iris.cube.Cube
-            Input cube to fix.
-
-        Returns
-        -------
-        iris.cube.Cube
-
-        """
-        rlat = cube.coord('grid_latitude').points
-        rlon = cube.coord('grid_longitude').points
-
-        # Transform grid latitude/longitude to array indices [0, 1, 2, ...]
-        rlat_to_idx = InterpolatedUnivariateSpline(rlat,
-                                                   np.arange(len(rlat)),
-                                                   k=1)
-        rlon_to_idx = InterpolatedUnivariateSpline(rlon,
-                                                   np.arange(len(rlon)),
-                                                   k=1)
-        rlat_idx_bnds = rlat_to_idx(cube.coord('grid_latitude').bounds)
-        rlon_idx_bnds = rlon_to_idx(cube.coord('grid_longitude').bounds)
-
-        # Calculate latitude/longitude vertices by interpolation
-        lat_vertices = []
-        lon_vertices = []
-        for (i, j) in [(0, 0), (0, 1), (1, 1), (1, 0)]:
-            (rlat_v, rlon_v) = np.meshgrid(rlat_idx_bnds[:, i],
-                                           rlon_idx_bnds[:, j],
-                                           indexing='ij')
-            lat_vertices.append(
-                map_coordinates(cube.coord('latitude').points,
-                                [rlat_v, rlon_v],
-                                mode='nearest'))
-            lon_vertices.append(
-                map_coordinates(cube.coord('longitude').points,
-                                [rlat_v, rlon_v],
-                                mode='wrap'))
-        lat_vertices = np.array(lat_vertices)
-        lon_vertices = np.array(lon_vertices)
-        lat_vertices = np.moveaxis(lat_vertices, 0, -1)
-        lon_vertices = np.moveaxis(lon_vertices, 0, -1)
-
-        # Copy vertices to cube
-        cube.coord('latitude').bounds = lat_vertices
-        cube.coord('longitude').bounds = lon_vertices
-
-        return cube
+Tos = OceanFixGrid

esmvalcore/cmor/_fixes/cmip5/bcc_csm1_1_m.py

@@ -1,9 +1,7 @@
 """Fixes for bcc-csm1-1-m."""
-from ..common import ClFixHybridPressureCoord
-from .bcc_csm1_1 import Tos as BaseTos
-
+from ..common import ClFixHybridPressureCoord, OceanFixGrid
 
 Cl = ClFixHybridPressureCoord
 
 
-Tos = BaseTos
+Tos = OceanFixGrid

esmvalcore/cmor/_fixes/cmip6/bcc_csm2_mr.py

@@ -1,7 +1,5 @@
 """Fixes for BCC-CSM2-MR model."""
-from ..cmip5.bcc_csm1_1 import Tos as BaseTos
-from ..common import ClFixHybridPressureCoord
-
+from ..common import ClFixHybridPressureCoord, OceanFixGrid
 
 Cl = ClFixHybridPressureCoord
 
@@ -12,57 +10,7 @@
 Clw = ClFixHybridPressureCoord
 
 
-class Tos(BaseTos):
-    """Fixes for tos."""
-
-    def fix_metadata(self, cubes):
-        """Rename ``var_name`` of 1D-``latitude`` and 1D-``longitude``.
-        Parameters
-        ----------
-        cubes : iris.cube.CubeList
-            Input cubes.
-        Returns
-        -------
-        iris.cube.CubeList
-        """
-        cube = self.get_cube_from_list(cubes)
-        lat_coord = cube.coord('latitude', dimensions=(1, ))
-        lon_coord = cube.coord('longitude', dimensions=(2, ))
-        lat_coord.standard_name = None
-        lat_coord.long_name = 'grid_latitude'
-        lat_coord.var_name = 'i'
-        lat_coord.units = '1'
-        lon_coord.standard_name = None
-        lon_coord.long_name = 'grid_longitude'
-        lon_coord.var_name = 'j'
-        lon_coord.units = '1'
-        lon_coord.circular = False
-        return cubes
-
+Tos = OceanFixGrid
 
-class Siconc(BaseTos):
-    """Fixes for siconc."""
 
-    def fix_metadata(self, cubes):
-        """Rename ``var_name`` of 1D-``latitude`` and 1D-``longitude``.
-        Parameters
-        ----------
-        cubes : iris.cube.CubeList
-            Input cubes.
-        Returns
-        -------
-        iris.cube.CubeList
-        """
-        cube = self.get_cube_from_list(cubes)
-        lat_coord = cube.coord('latitude', dimensions=(1, ))
-        lon_coord = cube.coord('longitude', dimensions=(2, ))
-        lat_coord.standard_name = None
-        lat_coord.long_name = 'grid_latitude'
-        lat_coord.var_name = 'i'
-        lat_coord.units = '1'
-        lon_coord.standard_name = None
-        lon_coord.long_name = 'grid_longitude'
-        lon_coord.var_name = 'j'
-        lon_coord.units = '1'
-        lon_coord.circular = False
-        return cubes
+Siconc = OceanFixGrid

esmvalcore/cmor/_fixes/cmip6/bcc_esm1.py

@@ -1,7 +1,5 @@
 """Fixes for BCC-ESM1 model."""
-from ..common import ClFixHybridPressureCoord
-from .bcc_csm2_mr import Tos as BaseTos
-
+from ..common import ClFixHybridPressureCoord, OceanFixGrid
 
 Cl = ClFixHybridPressureCoord
 
@@ -12,4 +10,4 @@
 Clw = ClFixHybridPressureCoord
 
 
-Tos = BaseTos
+Tos = OceanFixGrid

esmvalcore/cmor/_fixes/cmip6/fgoals_g3.py

@@ -2,7 +2,6 @@
 from ..cmip5.fgoals_g2 import Cl as BaseCl
 from ..common import OceanFixGrid
 
-
 Cl = BaseCl
 
 
@@ -12,7 +11,31 @@
 Clw = BaseCl
 
 
-Tos = OceanFixGrid
+class Tos(OceanFixGrid):
+    """Fixes for tos."""
+
+    def fix_metadata(self, cubes):
+        """Fix metadata.
+
+        FGOALS-g3 data contain latitude and longitude data set to >1e30 in some
+        places.
+
+        Parameters
+        ----------
+        cubes : iris.cube.CubeList
+            Input cubes.
+
+        Returns
+        -------
+        iris.cube.CubeList
+
+        """
+        cube = self.get_cube_from_list(cubes)
+        cube.coord('latitude').points[
+            cube.coord('latitude').points > 1000.0] = 0.0
+        cube.coord('longitude').points[
+            cube.coord('longitude').points > 1000.0] = 0.0
+        return super().fix_metadata(cubes)
 
 
-Siconc = OceanFixGrid
+Siconc = Tos

esmvalcore/cmor/_fixes/common.py

@@ -1,11 +1,15 @@
 """Common fixes used for multiple datasets."""
+import logging
+
 import iris
 import numpy as np
 from scipy.ndimage import map_coordinates
 
 from .fix import Fix
 from .shared import add_plev_from_altitude, fix_bounds
 
+logger = logging.getLogger(__name__)
+
 
 class ClFixHybridHeightCoord(Fix):
     """Fixes for ``cl`` regarding hybrid sigma height coordinates."""
@@ -108,49 +112,74 @@ def fix_metadata(self, cubes):
 class OceanFixGrid(Fix):
     """Fixes for tos, siconc in FGOALS-g3."""
 
-    def fix_data(self, cube):
-        """
-        Fix data.
-
-        Calculate missing latitude/longitude boundaries using interpolation.
-        Based on a similar fix for BCC-CSM2-MR.
+    def fix_metadata(self, cubes):
+        """Fix ``latitude`` and ``longitude`` (metadata and bounds).
 
         Parameters
         ----------
-        cube: iris.cube.Cube
-            Input cube to fix.
+        cubes : iris.cube.CubeList
+            Input cubes.
 
         Returns
         -------
-        iris.cube.Cube
+        iris.cube.CubeList
+
         """
-        rlat = cube.coord('grid_latitude').points
-        rlon = cube.coord('grid_longitude').points
-
-        # Guess coordinate bounds in rlat, rlon (following BCC-CSM2-MR-1).
-        rlat_idx_bnds = np.zeros((len(rlat), 2))
-        rlat_idx_bnds[:, 0] = np.arange(len(rlat)) - 0.5
-        rlat_idx_bnds[:, 1] = np.arange(len(rlat)) + 0.5
-        rlat_idx_bnds[0, 0] = 0.
-        rlat_idx_bnds[len(rlat) - 1, 1] = len(rlat)
-        rlon_idx_bnds = np.zeros((len(rlon), 2))
-        rlon_idx_bnds[:, 0] = np.arange(len(rlon)) - 0.5
-        rlon_idx_bnds[:, 1] = np.arange(len(rlon)) + 0.5
-
-        # Calculate latitude/longitude vertices by interpolation
+        cube = self.get_cube_from_list(cubes)
+        if cube.ndim != 3:
+            logger.warning(
+                "OceanFixGrid is designed to work on any data with an "
+                "irregular ocean grid, but it was only tested on 3D (time, "
+                "latitude, longitude) data so far; got %dD data", cube.ndim)
+
+        # Get dimensional coordinates. Note:
+        # - First dimension i -> X-direction (= longitude)
+        # - Second dimension j -> Y-direction (= latitude)
+        (j_dim, i_dim) = sorted(set(
+            cube.coord_dims(cube.coord('latitude', dim_coords=False)) +
+            cube.coord_dims(cube.coord('longitude', dim_coords=False))
+        ))
+        i_coord = cube.coord(dim_coords=True, dimensions=i_dim)
+        j_coord = cube.coord(dim_coords=True, dimensions=j_dim)
+
+        # Fix metadata of coordinate i
+        i_coord.var_name = 'i'
+        i_coord.standard_name = None
+        i_coord.long_name = 'cell index along first dimension'
+        i_coord.units = '1'
+        i_coord.circular = False
+
+        # Fix metadata of coordinate j
+        j_coord.var_name = 'j'
+        j_coord.standard_name = None
+        j_coord.long_name = 'cell index along second dimension'
+        j_coord.units = '1'
+
+        # Fix points and bounds of index coordinates i and j
+        for idx_coord in (i_coord, j_coord):
+            idx_coord.points = np.arange(len(idx_coord.points))
+            idx_coord.bounds = None
+            idx_coord.guess_bounds()
+
+        # Calculate latitude/longitude vertices by interpolation.
+        # Following the CF conventions (see
+        # cfconventions.org/cf-conventions/cf-conventions.html#cell-boundaries)
+        # we go counter-clockwise around the cells and construct a grid of
+        # index values which are in turn used to interpolate longitudes and
+        # latitudes in the midpoints between the cell centers.
         lat_vertices = []
         lon_vertices = []
-        for (i, j) in [(0, 0), (0, 1), (1, 1), (1, 0)]:
-            (rlat_v, rlon_v) = np.meshgrid(rlat_idx_bnds[:, i],
-                                           rlon_idx_bnds[:, j],
-                                           indexing='ij')
+        for (j, i) in [(0, 0), (0, 1), (1, 1), (1, 0)]:
+            (j_v, i_v) = np.meshgrid(j_coord.bounds[:, j],
+                                     i_coord.bounds[:, i],
+                                     indexing='ij')
             lat_vertices.append(
                 map_coordinates(cube.coord('latitude').points,
-                                [rlat_v, rlon_v],
+                                [j_v, i_v],
                                 mode='nearest'))
             lon_vertices.append(
                 map_coordinates(cube.coord('longitude').points,
-                                [rlat_v, rlon_v],
+                                [j_v, i_v],
                                 mode='wrap'))
         lat_vertices = np.array(lat_vertices)
         lon_vertices = np.array(lon_vertices)
@@ -160,39 +189,5 @@ def fix_data(self, cube):
         # Copy vertices to cube
         cube.coord('latitude').bounds = lat_vertices
         cube.coord('longitude').bounds = lon_vertices
-        return cube
-
-    def fix_metadata(self, cubes):
-        """
-        Rename ``var_name`` of 1D-``latitude`` and 1D-``longitude``.
-
-        Parameters
-        ----------
-        cubes : iris.cube.CubeList
-            Input cubes.
 
-        Returns
-        -------
-        iris.cube.CubeList
-        """
-        cube = self.get_cube_from_list(cubes)
-        lat_coord = cube.coord('cell index along second dimension',
-                               dimensions=(1, ))
-        lon_coord = cube.coord('cell index along first dimension',
-                               dimensions=(2, ))
-        lat_coord.standard_name = None
-        lat_coord.long_name = 'grid_latitude'
-        lat_coord.var_name = 'i'
-        lat_coord.units = '1'
-        lon_coord.standard_name = None
-        lon_coord.long_name = 'grid_longitude'
-        lon_coord.var_name = 'j'
-        lon_coord.units = '1'
-        lon_coord.circular = False
-        # FGOALS-g3 data contain latitude and longitude data set to
-        # >1e30 in some places. Set to 0. to avoid problem in check.py.
-        cube.coord('latitude').points[cube.coord('latitude').points > 1000.]\
-            = 0.
-        cube.coord('longitude').points[cube.coord('longitude').points > 1000.]\
-            = 0.
-        return cubes
+        return iris.cube.CubeList([cube])