Add general CORDEX coordinate boundary fix

esmvalcore/_recipe.py

@@ -337,6 +337,9 @@ def _add_fxvar_keys(fx_info, variable):
     # add special ensemble for CMIP5 only
     if fx_variable['project'] == 'CMIP5':
         fx_variable['ensemble'] = 'r0i0p0'
+    elif fx_variable['project'] in ['CORDEX']:
+        fx_variable['ensemble'] = [fx_variable['ensemble'], 'r0i0p0']
+
 
     # add missing cmor info
     _add_cmor_info(fx_variable, override=True)

esmvalcore/cmor/_fixes/cordex/__init__.py

@@ -0,0 +1 @@
+"""Fixes for CORDEX data."""

esmvalcore/cmor/_fixes/cordex/project.py

@@ -0,0 +1,251 @@
+"""Fixes for CORDEX project."""
+import logging
+import numpy as np
+import iris
+
+from esmvalcore.cmor.fix import Fix
+from esmvalcore.preprocessor._shared import guess_bounds
+
+logger = logging.getLogger(__name__)
+
+
+class AllVars(Fix):
+    """Generic fixes for the CORDEX project."""
+
+    def fix_metadata(self, cubes):
+        """Fix metatdata.
+
+        Calculate missing latitude/longitude boundaries using interpolation.
+
+        Parameters
+        ----------
+        cubes: iris.cube.CubeList
+            List of cubes to fix
+
+        Returns
+        -------
+        iris.cube.CubeList
+            Fixed cubes
+
+        """
+        cube = self.get_cube_from_list(cubes)
+
+        coord_names = [coord.standard_name for coord in cube.coords()]
+        x_name = cube.coord(axis='x', dim_coords=True).standard_name
+        y_name = cube.coord(axis='y', dim_coords=True).standard_name
+        if x_name != 'longitude' and y_name != 'latitude':
+            # Fix some issues with Lambert conformal coordinates
+            self.fix_coordinate_system(cube)
+
+            # get the grid bounds
+            cube = guess_bounds(cube, [x_name, y_name])
+
+            # get the lat lon the for global coord system (if missing)
+            if 'latitude' not in coord_names or \
+               'longitude' not in coord_names:
+                cube = self.get_lat_lon_coordinates(cube)
+            # get the lat lon bounds for the global coord system (if missing)
+            if not cube.coord('latitude').has_bounds() and \
+               not cube.coord('longitude').has_bounds():
+                cube = self.get_lat_lon_bounds(cube)
+
+            cube = self.check_grid_differences(cube)
+
+        return cubes
+
+    @staticmethod
+    def check_grid_differences(cube):
+        """
+        Derive lat and lon coordinates from grid coordinates.
+
+        It also warns about the maximum differences
+        """
+        x_coord = cube.coord(axis='x', dim_coords=True)
+        y_coord = cube.coord(axis='y', dim_coords=True)
+        glon, glat = np.meshgrid(x_coord.points, y_coord.points)
+
+        global_coord_sys = iris.coord_systems.GeogCS(
+            iris.fileformats.pp.EARTH_RADIUS)
+        grid_coord_sys = x_coord.coord_system
+
+        if isinstance(grid_coord_sys, iris.coord_systems.RotatedGeogCS):
+            gr_type = "rotated"
+            glon, glat = np.meshgrid(x_coord.points, y_coord.points)
+            lons, lats = iris.analysis.cartography.unrotate_pole(
+                glon, glat,
+                grid_coord_sys.grid_north_pole_longitude,
+                grid_coord_sys.grid_north_pole_latitude
+            )
+        elif isinstance(grid_coord_sys, iris.coord_systems.LambertConformal):
+            gr_type = "lcc"
+            ccrs_global = global_coord_sys.as_cartopy_crs()
+            xyz = ccrs_global.transform_points(grid_coord_sys.as_cartopy_crs(),
+                                               glon, glat)
+            lons = xyz[:, :, 0]
+            lats = xyz[:, :, 1]
+        else:
+            emsg = 'Unknown coordinate system, got {!r}.'
+            raise NotImplementedError(
+                emsg.format(type(x_coord.grid_coord_sys)))
+
+        dif = np.sqrt(np.square(lats - cube.coord('latitude').points) +
+                      np.square(lons - cube.coord('longitude').points))
+        if dif.max() != 0:
+            logger.warning("There are diffs. betw. %s and lat-lon.", gr_type)
+            logger.warning("Max diff: %s ; Min diff: %s", dif.max(), dif.min())
+            cube.attributes.update({'grid_max_spacial_diff': dif.max()})
+            cube.attributes.update({'grid_min_spacial_diff': dif.min()})
+        return cube
+
+    @staticmethod
+    def get_lat_lon_bounds(cube):
+        """
+        Derive lat and lon bounds from grid coordinates.
+
+        CMOR standard for 2-D coordinate bounds is indexing the four vertices
+        as follows: 0=(j-1,i-1), 1=(j-1,i+1), 2=(j+1,i+1), 3=(j+1,i-1).
+        """
+        x_coord = cube.coord(axis='x', dim_coords=True)
+        y_coord = cube.coord(axis='y', dim_coords=True)
+        glon, glat = np.meshgrid(x_coord.bounds, y_coord.bounds)
+
+        global_coord_sys = iris.coord_systems.GeogCS(
+            iris.fileformats.pp.EARTH_RADIUS)
+        grid_coord_sys = x_coord.coord_system
+
+        if isinstance(grid_coord_sys, iris.coord_systems.RotatedGeogCS):
+            lon_bnds, lat_bnds = iris.analysis.cartography.unrotate_pole(
+                glon, glat,
+                grid_coord_sys.grid_north_pole_longitude,
+                grid_coord_sys.grid_north_pole_latitude
+            )
+        elif isinstance(grid_coord_sys, iris.coord_systems.LambertConformal):
+            ccrs_global = global_coord_sys.as_cartopy_crs()
+            xyz = ccrs_global.transform_points(
+                grid_coord_sys.as_cartopy_crs(), glon, glat)
+            lon_bnds = xyz[:, :, 0]
+            lat_bnds = xyz[:, :, 1]
+        else:
+            emsg = 'Unknown coordinate system, got {!r}.'
+            raise NotImplementedError(
+                emsg.format(type(x_coord.grid_coord_sys)))
+
+        if not cube.coord('latitude').has_bounds():
+            lat_bnds = [lat_bnds[::2, ::2], lat_bnds[::2, 1::2],
+                        lat_bnds[1::2, 1::2], lat_bnds[1::2, ::2]]
+            lat_bnds = np.array(lat_bnds)
+            lat_bnds = np.moveaxis(lat_bnds, 0, -1)
+            cube.coord('latitude').bounds = lat_bnds
+
+        if not cube.coord('longitude').has_bounds():
+            lon_bnds = [lon_bnds[::2, ::2], lon_bnds[::2, 1::2],
+                        lon_bnds[1::2, 1::2], lon_bnds[1::2, ::2]]
+            lon_bnds = np.array(lon_bnds)
+            lon_bnds = np.moveaxis(lon_bnds, 0, -1)
+            cube.coord('longitude').bounds = lon_bnds
+
+        return cube
+
+    @staticmethod
+    def get_lat_lon_coordinates(cube):
+        """Derive lat and lon coordinates from grid coordinates."""
+        coord_names = [coord.standard_name for coord in cube.coords()]
+        x_coord = cube.coord(axis='x', dim_coords=True)
+        y_coord = cube.coord(axis='y', dim_coords=True)
+        glon, glat = np.meshgrid(x_coord.points, y_coord.points)
+
+        global_coord_sys = iris.coord_systems.GeogCS(
+            iris.fileformats.pp.EARTH_RADIUS)
+        grid_coord_sys = x_coord.coord_system
+
+        if isinstance(grid_coord_sys, iris.coord_systems.RotatedGeogCS):
+            glon, glat = np.meshgrid(x_coord.points, y_coord.points)
+            lons, lats = iris.analysis.cartography.unrotate_pole(
+                glon, glat,
+                grid_coord_sys.grid_north_pole_longitude,
+                grid_coord_sys.grid_north_pole_latitude
+            )
+        elif isinstance(grid_coord_sys, iris.coord_systems.LambertConformal):
+            ccrs_global = global_coord_sys.as_cartopy_crs()
+            xyz = ccrs_global.transform_points(
+                grid_coord_sys.as_cartopy_crs(), glon, glat)
+            lons = xyz[:, :, 0]
+            lats = xyz[:, :, 1]
+        else:
+            emsg = 'Unknown coordinate system, got {!r}.'
+            raise NotImplementedError(
+                emsg.format(type(x_coord.grid_coord_sys)))
+
+        if 'latitude' not in coord_names:
+            lat = iris.coords.AuxCoord(lats,
+                                       var_name='lat',
+                                       standard_name='latitude',
+                                       units='degrees',
+                                       coord_system=global_coord_sys)
+            AllVars._add_coordinate(lat, cube)
+
+        if 'longitude' not in coord_names:
+            lon = iris.coords.AuxCoord(lons,
+                                       var_name='lon',
+                                       standard_name='longitude',
+                                       units='degrees',
+                                       coord_system=global_coord_sys)
+            AllVars._add_coordinate(lon, cube)
+
+        return cube
+
+    @staticmethod
+    def _add_coordinate(lat, cube):
+        if lat.shape[0] == lat.shape[1]:
+            data_dims = np.where(np.array(cube.shape) == lat.shape[0])[0]
+        else:
+            data_dims = [cube.shape.index(lat.shape[0]),
+                         cube.shape.index(lat.shape[1])]
+
+        cube.add_aux_coord(lat, data_dims=data_dims)
+
+    @staticmethod
+    def fix_coordinate_system(cube):
+        """Fix LambertConformal."""
+        x_coord = cube.coord(axis='x', dim_coords=True)
+        y_coord = cube.coord(axis='y', dim_coords=True)
+
+        global_coord_sys = iris.coord_systems.GeogCS(
+            iris.fileformats.pp.EARTH_RADIUS)
+        grid_coord_sys = x_coord.coord_system
+
+        if isinstance(grid_coord_sys, iris.coord_systems.LambertConformal):
+            # there are some badly defined coordinate systems:
+            if grid_coord_sys.ellipsoid:
+                if grid_coord_sys.ellipsoid.semi_minor_axis == float('-inf'):
+                    grid_coord_sys.ellipsoid.semi_minor_axis = None
+            else:
+                grid_coord_sys.ellipsoid = global_coord_sys
+
+            # badly defined units
+            if x_coord.units == 'km':
+                x_coord.convert_units('meter')
+            if y_coord.units == 'km':
+                y_coord.convert_units('meter')
+
+            # and badly defined offsets
+            if x_coord[0].points[0] == 0:
+                x_offset = x_coord.points.mean()
+                grid_coord_sys.false_easting = x_offset
+            if y_coord[0].points[0] == 0:
+                y_offset = y_coord.points.mean()
+                grid_coord_sys.false_northing = y_offset
+
+            # This is for
+            # - {dataset: ICTP-RegCM4-3, driver: MOHC-HadGEM2-ES}
+            # LambertConformal of the netcdf gives a false_easting though
+            # x_offset is 0, setting false_easting to zero gives much smaller
+            # differences to latitude longitude in the netcdf
+            if x_coord.points.mean() == 0 and \
+               grid_coord_sys.false_easting is not None:
+                grid_coord_sys.false_easting = 0
+            if y_coord.points.mean() == 0 and \
+               grid_coord_sys.false_northing is not None:
+                grid_coord_sys.false_northing = 0
+
+        return cube

esmvalcore/cmor/_fixes/fix.py

@@ -1,7 +1,7 @@
 """Contains the base class for dataset fixes."""
 import importlib
-import os
 import inspect
+import os
 
 from ..table import CMOR_TABLES
 
@@ -157,19 +157,20 @@ def get_fixes(project, dataset, mip, short_name, extra_facets=None):
             extra_facets = {}
 
         fixes = []
-        try:
-            fixes_module = importlib.import_module(
-                'esmvalcore.cmor._fixes.{0}.{1}'.format(project, dataset))
-
-            classes = inspect.getmembers(fixes_module, inspect.isclass)
-            classes = dict((name.lower(), value) for name, value in classes)
-            for fix_name in (short_name, mip.lower(), 'allvars'):
-                try:
-                    fixes.append(classes[fix_name](vardef, extra_facets))
-                except KeyError:
-                    pass
-        except ImportError:
-            pass
+        for package in ('project', dataset):
+            try:
+                fixes_module = importlib.import_module(
+                    f'esmvalcore.cmor._fixes.{project}.{package}')
+
+                classes = inspect.getmembers(fixes_module, inspect.isclass)
+                classes = dict((name.lower(), val) for name, val in classes)
+                for fix_name in (short_name, mip.lower(), 'allvars'):
+                    try:
+                        fixes.append(classes[fix_name](vardef))
+                    except KeyError:
+                        pass
+            except ImportError:
+                pass
         return fixes
 
     @staticmethod

tests/integration/cmor/_fixes/cordex/test_generic_cordex.py

@@ -0,0 +1,55 @@
+"""Tests for the fixes of CORDEX."""
+import iris
+import iris.coord_systems
+from iris.coords import AuxCoord, DimCoord
+
+from esmvalcore.cmor._fixes.cordex.project import AllVars
+from esmvalcore.cmor.fix import Fix
+
+
+def test_get_allvars_fix():
+    fix = Fix.get_fixes('CORDEX', 'any_dataset', 'mip', 'tas')
+    assert fix == [AllVars(None)]
+
+
+def test_allvars():
+    cube = iris.cube.Cube(
+        [[1.0, 1.0], [1.0, 1.0]],
+        var_name='co2',
+        standard_name='mole_fraction_of_carbon_dioxide_in_air',
+        units='mol mol-1',
+    )
+    # Rotated grids
+    grid_latitude = DimCoord(, units=degress)
+    grid_longitude = DimCoord(, units=degress)
+    latitude = AuxCoord()
+    longitude = AuxCoord()
+    coord_sys_rotated = iris.coord_systems.RotatedGeogCS(39.25, -162.0)
+    grid_lat_11 = iris.coords.DimCoord(np.linspace(-23.375, 21.835, 412),
+                                    var_name='rlat',
+                                    standard_name='grid_latitude',
+                                    long_name='latitude in rotated-pole grid',
+                                    units='degrees',
+                                    coord_system=coord_sys_rotated)
+    grid_lon_11 = iris.coords.DimCoord(np.linspace(-28.375, 18.155, 424),
+                                    var_name='rlon',
+                                    standard_name='grid_longitude',
+                                    long_name='longitude in rotated-pole grid',
+                                    units='degrees',
+                                    coord_system=coord_sys_rotated)
+
+
+def test_allvars_lambert_conformal():
+    cube = iris.cube.Cube(
+        [[1.0, 1.0], [1.0, 1.0]],
+        var_name='co2',
+        standard_name='mole_fraction_of_carbon_dioxide_in_air',
+        units='mol mol-1',
+    )
+    # Rotated grids
+    projection_x_coordinate = DimCoord([-1000, 1000],units=m, coord_system=iris.coord_systems.LambertConformal)
+
+
+    LambertConformal(central_lat=48.0, central_lon=9.75, false_easting=-6000.0, false_northing=-6000.0, secant_latitudes=(30.0, 65.0), ellipsoid=GeogCS(semi_major_axis=6371229.0, semi_minor_axis=-inf))
+
+    grid_longitude = DimCoord()