diff --git a/esmvalcore/_recipe.py b/esmvalcore/_recipe.py index 8d361efc62..c66f7a1f62 100644 --- a/esmvalcore/_recipe.py +++ b/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) diff --git a/esmvalcore/cmor/_fixes/cordex/__init__.py b/esmvalcore/cmor/_fixes/cordex/__init__.py new file mode 100644 index 0000000000..093969370f --- /dev/null +++ b/esmvalcore/cmor/_fixes/cordex/__init__.py @@ -0,0 +1 @@ +"""Fixes for CORDEX data.""" diff --git a/esmvalcore/cmor/_fixes/cordex/project.py b/esmvalcore/cmor/_fixes/cordex/project.py new file mode 100644 index 0000000000..61e5c96a53 --- /dev/null +++ b/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 diff --git a/esmvalcore/cmor/_fixes/fix.py b/esmvalcore/cmor/_fixes/fix.py index be94c4a774..bfbf30fb08 100644 --- a/esmvalcore/cmor/_fixes/fix.py +++ b/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 diff --git a/tests/integration/cmor/_fixes/cordex/test_generic_cordex.py b/tests/integration/cmor/_fixes/cordex/test_generic_cordex.py new file mode 100644 index 0000000000..63a35dd177 --- /dev/null +++ b/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()