From 358b7048ba84ba9c598600957ea31a92646fb19b Mon Sep 17 00:00:00 2001 From: aperez Date: Tue, 23 Jul 2019 11:25:37 +0200 Subject: [PATCH 01/11] Add general CORDEX coordinate boundary fix --- .../CORDEX/Cordex_boundaries_guessing.py | 60 +++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 esmvalcore/cmor/_fixes/CORDEX/Cordex_boundaries_guessing.py diff --git a/esmvalcore/cmor/_fixes/CORDEX/Cordex_boundaries_guessing.py b/esmvalcore/cmor/_fixes/CORDEX/Cordex_boundaries_guessing.py new file mode 100644 index 0000000000..1f7f2a31e5 --- /dev/null +++ b/esmvalcore/cmor/_fixes/CORDEX/Cordex_boundaries_guessing.py @@ -0,0 +1,60 @@ +"""Fixes for bcc-csm1-1.""" +import numpy as np +from scipy.interpolate import InterpolatedUnivariateSpline +from scipy.ndimage import map_coordinates + +from ..fix import Fix + + +class CordexFix(Fix): + """Fixes for tos.""" + + def fix_data(self, cube): + """Fix data. + + Calculate missing latitude/longitude boundaries using interpolation. + + Parameters + ---------- + cube: iris.cube.Cube + + 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')._guess_bounds()) + rlon_idx_bnds = rlon_to_idx(cube.coord('grid_longitude')._guess_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 From a5693947f3219a6708bf37faa6d487f6d8cb796f Mon Sep 17 00:00:00 2001 From: Javier Vegas-Regidor Date: Tue, 5 May 2020 08:06:47 +0200 Subject: [PATCH 02/11] Add generic cordex fix --- .../CORDEX/Cordex_boundaries_guessing.py | 60 ---------------- esmvalcore/cmor/_fixes/cordex/__init__.py | 1 + esmvalcore/cmor/_fixes/cordex/project.py | 70 +++++++++++++++++++ esmvalcore/cmor/_fixes/fix.py | 27 +++---- .../cmor/_fixes/cordex/test_generic_cordex.py | 9 +++ 5 files changed, 94 insertions(+), 73 deletions(-) delete mode 100644 esmvalcore/cmor/_fixes/CORDEX/Cordex_boundaries_guessing.py create mode 100644 esmvalcore/cmor/_fixes/cordex/__init__.py create mode 100644 esmvalcore/cmor/_fixes/cordex/project.py create mode 100644 tests/integration/cmor/_fixes/cordex/test_generic_cordex.py diff --git a/esmvalcore/cmor/_fixes/CORDEX/Cordex_boundaries_guessing.py b/esmvalcore/cmor/_fixes/CORDEX/Cordex_boundaries_guessing.py deleted file mode 100644 index 1f7f2a31e5..0000000000 --- a/esmvalcore/cmor/_fixes/CORDEX/Cordex_boundaries_guessing.py +++ /dev/null @@ -1,60 +0,0 @@ -"""Fixes for bcc-csm1-1.""" -import numpy as np -from scipy.interpolate import InterpolatedUnivariateSpline -from scipy.ndimage import map_coordinates - -from ..fix import Fix - - -class CordexFix(Fix): - """Fixes for tos.""" - - def fix_data(self, cube): - """Fix data. - - Calculate missing latitude/longitude boundaries using interpolation. - - Parameters - ---------- - cube: iris.cube.Cube - - 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')._guess_bounds()) - rlon_idx_bnds = rlon_to_idx(cube.coord('grid_longitude')._guess_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 diff --git a/esmvalcore/cmor/_fixes/cordex/__init__.py b/esmvalcore/cmor/_fixes/cordex/__init__.py new file mode 100644 index 0000000000..290b07295f --- /dev/null +++ b/esmvalcore/cmor/_fixes/cordex/__init__.py @@ -0,0 +1 @@ +"""Fixes for CORDEX data""" \ No newline at end of file diff --git a/esmvalcore/cmor/_fixes/cordex/project.py b/esmvalcore/cmor/_fixes/cordex/project.py new file mode 100644 index 0000000000..de3d919d29 --- /dev/null +++ b/esmvalcore/cmor/_fixes/cordex/project.py @@ -0,0 +1,70 @@ +"""Fixes for CORDEX project""" +import numpy as np +from scipy.interpolate import InterpolatedUnivariateSpline +from scipy.ndimage import map_coordinates + +from esmvalcore.cmor.fix import Fix + + +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 + + Returns + ------- + iris.cube.CubeList + + """ + for cube in cubes: + if not cube.coords('latitude'): + continue + if not cube.coords('longitude'): + continue + if cube.coord('latitude').has_bounds() and \ + cube.coord('longitude').has_bounds(): + continue + 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').guess_bounds()) + rlon_idx_bnds = rlon_to_idx( + cube.coord('grid_longitude').guess_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 cubes diff --git a/esmvalcore/cmor/_fixes/fix.py b/esmvalcore/cmor/_fixes/fix.py index 12aee7b9a0..c53e25e8ff 100644 --- a/esmvalcore/cmor/_fixes/fix.py +++ b/esmvalcore/cmor/_fixes/fix.py @@ -156,19 +156,20 @@ def get_fixes(project, dataset, mip, short_name): short_name = short_name.replace('-', '_').lower() 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, 'allvars'): - try: - fixes.append(classes[fix_name](vardef)) - 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(), value) for name, value in classes) + for fix_name in (short_name, '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..89d5ee0949 --- /dev/null +++ b/tests/integration/cmor/_fixes/cordex/test_generic_cordex.py @@ -0,0 +1,9 @@ +"""Tests for the fixes of CORDEX.""" + +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', 'tas') + assert fix == [AllVars()] From 2602bfba810620ec74056069b09e6cc556d84ecb Mon Sep 17 00:00:00 2001 From: Javier Vegas-Regidor Date: Tue, 5 May 2020 08:53:36 +0200 Subject: [PATCH 03/11] Fix tests --- esmvalcore/cmor/_fixes/cordex/__init__.py | 2 +- tests/integration/cmor/_fixes/cordex/test_generic_cordex.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/esmvalcore/cmor/_fixes/cordex/__init__.py b/esmvalcore/cmor/_fixes/cordex/__init__.py index 290b07295f..02d8843142 100644 --- a/esmvalcore/cmor/_fixes/cordex/__init__.py +++ b/esmvalcore/cmor/_fixes/cordex/__init__.py @@ -1 +1 @@ -"""Fixes for CORDEX data""" \ No newline at end of file +"""Fixes for CORDEX data""" diff --git a/tests/integration/cmor/_fixes/cordex/test_generic_cordex.py b/tests/integration/cmor/_fixes/cordex/test_generic_cordex.py index 89d5ee0949..ef3954518b 100644 --- a/tests/integration/cmor/_fixes/cordex/test_generic_cordex.py +++ b/tests/integration/cmor/_fixes/cordex/test_generic_cordex.py @@ -5,5 +5,5 @@ def test_get_allvars_fix(): - fix = Fix.get_fixes('CORDEX', 'any_dataset', 'tas') - assert fix == [AllVars()] + fix = Fix.get_fixes('CORDEX', 'any_dataset', 'mip', 'tas') + assert fix == [AllVars(None)] From 07c9e8b6b720c0b333063b3daf7d153a8a34bc8e Mon Sep 17 00:00:00 2001 From: Javier Vegas-Regidor Date: Tue, 5 May 2020 13:05:25 +0200 Subject: [PATCH 04/11] Fix cordex --- esmvalcore/cmor/_fixes/cordex/project.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/esmvalcore/cmor/_fixes/cordex/project.py b/esmvalcore/cmor/_fixes/cordex/project.py index de3d919d29..2cf1a2147b 100644 --- a/esmvalcore/cmor/_fixes/cordex/project.py +++ b/esmvalcore/cmor/_fixes/cordex/project.py @@ -40,9 +40,9 @@ def fix_metadata(self, cubes): rlon_to_idx = InterpolatedUnivariateSpline( rlon, np.arange(len(rlon)), k=1) rlat_idx_bnds = rlat_to_idx( - cube.coord('grid_latitude').guess_bounds()) + cube.coord('grid_latitude')._guess_bounds()) rlon_idx_bnds = rlon_to_idx( - cube.coord('grid_longitude').guess_bounds()) + cube.coord('grid_longitude')._guess_bounds()) # Calculate latitude/longitude vertices by interpolation lat_vertices = [] From f11582dca7e92a3831d02ae9f65166ecd0806dbc Mon Sep 17 00:00:00 2001 From: Javier Vegas-Regidor Date: Wed, 6 May 2020 17:50:43 +0200 Subject: [PATCH 05/11] Avoid using private functions --- esmvalcore/cmor/_fixes/cordex/project.py | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/esmvalcore/cmor/_fixes/cordex/project.py b/esmvalcore/cmor/_fixes/cordex/project.py index 2cf1a2147b..3e9f53ed89 100644 --- a/esmvalcore/cmor/_fixes/cordex/project.py +++ b/esmvalcore/cmor/_fixes/cordex/project.py @@ -31,18 +31,11 @@ def fix_metadata(self, cubes): if cube.coord('latitude').has_bounds() and \ cube.coord('longitude').has_bounds(): continue - 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')._guess_bounds()) - rlon_idx_bnds = rlon_to_idx( - cube.coord('grid_longitude')._guess_bounds()) + rlat_idx_bnds = self._get_points_and_bounds( + cube.coord('grid_latitude')) + rlon_idx_bnds = self._get_points_and_bounds( + cube.coord('grid_longitude')) # Calculate latitude/longitude vertices by interpolation lat_vertices = [] @@ -68,3 +61,11 @@ def fix_metadata(self, cubes): cube.coord('longitude').bounds = lon_vertices return cubes + + def _get_points_and_bounds(self, coord): + coord = coord.copy() + coord.guess_bounds() + points_to_idx = InterpolatedUnivariateSpline( + coord.points, np.arange(len(coord.points)), k=1) + bounds = points_to_idx(coord.bounds) + return bounds From 9b021a7da777d9bda6285129d229816301c1b52c Mon Sep 17 00:00:00 2001 From: Javier Vegas-Regidor Date: Thu, 1 Oct 2020 16:41:28 +0200 Subject: [PATCH 06/11] Update project.py --- esmvalcore/cmor/_fixes/cordex/project.py | 259 +++++++++++++++++++---- esmvalcore/cmor/_fixes/fix.py | 2 +- 2 files changed, 215 insertions(+), 46 deletions(-) diff --git a/esmvalcore/cmor/_fixes/cordex/project.py b/esmvalcore/cmor/_fixes/cordex/project.py index 3e9f53ed89..95fd9529e9 100644 --- a/esmvalcore/cmor/_fixes/cordex/project.py +++ b/esmvalcore/cmor/_fixes/cordex/project.py @@ -1,9 +1,12 @@ """Fixes for CORDEX project""" +import logging import numpy as np -from scipy.interpolate import InterpolatedUnivariateSpline -from scipy.ndimage import map_coordinates +import iris from esmvalcore.cmor.fix import Fix +from esmvalcore.preprocessor._shared import guess_bounds + +logger = logging.getLogger(__name__) class AllVars(Fix): @@ -23,49 +26,215 @@ def fix_metadata(self, cubes): iris.cube.CubeList """ - for cube in cubes: - if not cube.coords('latitude'): - continue - if not cube.coords('longitude'): - continue - if cube.coord('latitude').has_bounds() and \ - cube.coord('longitude').has_bounds(): - continue - - rlat_idx_bnds = self._get_points_and_bounds( - cube.coord('grid_latitude')) - rlon_idx_bnds = self._get_points_and_bounds( - cube.coord('grid_longitude')) - - # 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 + 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 - def _get_points_and_bounds(self, coord): - coord = coord.copy() - coord.guess_bounds() - points_to_idx = InterpolatedUnivariateSpline( - coord.points, np.arange(len(coord.points)), k=1) - bounds = points_to_idx(coord.bounds) - return bounds + def check_grid_differences(self, cube): + """Derive lat and lon coordinates from grid coordinates. + And warn 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(f"There are diffs. betw. {gr_type} and lat-lon.") + logger.warning(f"Max diff: {dif.max()} ; Min diff: {dif.min()}") + cube.attributes.update({'grid_max_spacial_diff': dif.max()}) + cube.attributes.update({'grid_min_spacial_diff': dif.min()}) + return cube + + def get_lat_lon_bounds(self, cube): + """ + Derive lat and lon bounds from grid coordinates. + CMOR standard for 2-D coordinate counds 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 + + def get_lat_lon_coordinates(self, 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) + 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) + + if 'longitude' not in coord_names: + lon = iris.coords.AuxCoord(lons, + var_name='lon', + standard_name='longitude', + units='degrees', + coord_system=global_coord_sys) + if lon.shape[0] == lon.shape[1]: + data_dims = np.where(np.array(cube.shape) == lon.shape[0])[0] + else: + data_dims = [cube.shape.index(lon.shape[0]), + cube.shape.index(lon.shape[1])] + + cube.add_aux_coord(lon, data_dims=data_dims) + + return cube + + def fix_coordinate_system(self, 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 c53e25e8ff..e4b891bed5 100644 --- a/esmvalcore/cmor/_fixes/fix.py +++ b/esmvalcore/cmor/_fixes/fix.py @@ -162,7 +162,7 @@ def get_fixes(project, dataset, mip, short_name): f'esmvalcore.cmor._fixes.{project}.{package}') classes = inspect.getmembers(fixes_module, inspect.isclass) - classes = dict((name.lower(), value) for name, value in classes) + classes = dict((name.lower(), val) for name, val in classes) for fix_name in (short_name, 'allvars'): try: fixes.append(classes[fix_name](vardef)) From 93d69ed5ca7b08e5048338f0feb0301c270cf67d Mon Sep 17 00:00:00 2001 From: Javier Vegas-Regidor Date: Thu, 1 Oct 2020 18:10:46 +0200 Subject: [PATCH 07/11] Fix docstring warnings --- esmvalcore/cmor/_fixes/cordex/project.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/esmvalcore/cmor/_fixes/cordex/project.py b/esmvalcore/cmor/_fixes/cordex/project.py index 95fd9529e9..6fb61f4fe5 100644 --- a/esmvalcore/cmor/_fixes/cordex/project.py +++ b/esmvalcore/cmor/_fixes/cordex/project.py @@ -1,4 +1,4 @@ -"""Fixes for CORDEX project""" +"""Fixes for CORDEX project.""" import logging import numpy as np import iris @@ -10,7 +10,7 @@ class AllVars(Fix): - """Generic fixes for the CORDEX project""" + """Generic fixes for the CORDEX project.""" def fix_metadata(self, cubes): """Fix metatdata. @@ -52,8 +52,11 @@ def fix_metadata(self, cubes): return cubes def check_grid_differences(self, cube): - """Derive lat and lon coordinates from grid coordinates. - And warn about the maximum differences""" + """ + 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) @@ -93,7 +96,8 @@ def check_grid_differences(self, cube): def get_lat_lon_bounds(self, cube): """ Derive lat and lon bounds from grid coordinates. - CMOR standard for 2-D coordinate counds is indexing the four vertices + + 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) From a1e362ce751170a52e7b35a67a0f1b366d6d6f6f Mon Sep 17 00:00:00 2001 From: Javier Vegas-Regidor Date: Thu, 1 Oct 2020 18:17:55 +0200 Subject: [PATCH 08/11] Fix a bunch of codacy issues --- esmvalcore/cmor/_fixes/cordex/project.py | 45 ++++++++++++++---------- 1 file changed, 26 insertions(+), 19 deletions(-) diff --git a/esmvalcore/cmor/_fixes/cordex/project.py b/esmvalcore/cmor/_fixes/cordex/project.py index 6fb61f4fe5..32d1d8d30a 100644 --- a/esmvalcore/cmor/_fixes/cordex/project.py +++ b/esmvalcore/cmor/_fixes/cordex/project.py @@ -51,7 +51,8 @@ def fix_metadata(self, cubes): return cubes - def check_grid_differences(self, cube): + @staticmethod + def check_grid_differences(cube): """ Derive lat and lon coordinates from grid coordinates. @@ -69,9 +70,10 @@ def check_grid_differences(self, cube): 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) + 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() @@ -87,13 +89,14 @@ def check_grid_differences(self, cube): dif = np.sqrt(np.square(lats - cube.coord('latitude').points) + np.square(lons - cube.coord('longitude').points)) if dif.max() != 0: - logger.warning(f"There are diffs. betw. {gr_type} and lat-lon.") - logger.warning(f"Max diff: {dif.max()} ; Min diff: {dif.min()}") + 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 - def get_lat_lon_bounds(self, cube): + @staticmethod + def get_lat_lon_bounds(cube): """ Derive lat and lon bounds from grid coordinates. @@ -110,13 +113,14 @@ def get_lat_lon_bounds(self, cube): 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) + 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) + xyz = ccrs_global.transform_points( + grid_coord_sys.as_cartopy_crs(), glon, glat) lon_bnds = xyz[:, :, 0] lat_bnds = xyz[:, :, 1] else: @@ -140,7 +144,8 @@ def get_lat_lon_bounds(self, cube): return cube - def get_lat_lon_coordinates(self, 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) @@ -154,13 +159,14 @@ def get_lat_lon_coordinates(self, cube): 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) + 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) + xyz = ccrs_global.transform_points( + grid_coord_sys.as_cartopy_crs(), glon, glat) lons = xyz[:, :, 0] lats = xyz[:, :, 1] else: @@ -198,7 +204,8 @@ def get_lat_lon_coordinates(self, cube): return cube - def fix_coordinate_system(self, cube): + @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) From 4492581bd5d083431874dd4948fc932e5b658418 Mon Sep 17 00:00:00 2001 From: Javier Vegas-Regidor Date: Thu, 1 Oct 2020 18:18:50 +0200 Subject: [PATCH 09/11] Fix doc issue --- esmvalcore/cmor/_fixes/cordex/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/esmvalcore/cmor/_fixes/cordex/__init__.py b/esmvalcore/cmor/_fixes/cordex/__init__.py index 02d8843142..093969370f 100644 --- a/esmvalcore/cmor/_fixes/cordex/__init__.py +++ b/esmvalcore/cmor/_fixes/cordex/__init__.py @@ -1 +1 @@ -"""Fixes for CORDEX data""" +"""Fixes for CORDEX data.""" From 13738b101f5b6e928118ae11d459858fefc8e790 Mon Sep 17 00:00:00 2001 From: Javier Vegas-Regidor Date: Fri, 2 Oct 2020 10:25:06 +0200 Subject: [PATCH 10/11] Avoid duplication --- esmvalcore/cmor/_fixes/cordex/project.py | 28 ++++++++++++------------ 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/esmvalcore/cmor/_fixes/cordex/project.py b/esmvalcore/cmor/_fixes/cordex/project.py index 32d1d8d30a..61e5c96a53 100644 --- a/esmvalcore/cmor/_fixes/cordex/project.py +++ b/esmvalcore/cmor/_fixes/cordex/project.py @@ -20,10 +20,12 @@ def fix_metadata(self, cubes): Parameters ---------- cubes: iris.cube.CubeList + List of cubes to fix Returns ------- iris.cube.CubeList + Fixed cubes """ cube = self.get_cube_from_list(cubes) @@ -180,13 +182,7 @@ def get_lat_lon_coordinates(cube): standard_name='latitude', units='degrees', coord_system=global_coord_sys) - 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) + AllVars._add_coordinate(lat, cube) if 'longitude' not in coord_names: lon = iris.coords.AuxCoord(lons, @@ -194,16 +190,20 @@ def get_lat_lon_coordinates(cube): standard_name='longitude', units='degrees', coord_system=global_coord_sys) - if lon.shape[0] == lon.shape[1]: - data_dims = np.where(np.array(cube.shape) == lon.shape[0])[0] - else: - data_dims = [cube.shape.index(lon.shape[0]), - cube.shape.index(lon.shape[1])] - - cube.add_aux_coord(lon, data_dims=data_dims) + 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.""" From 2ed7ca3fc0d47fbbe6161a3842369db3fd4a7824 Mon Sep 17 00:00:00 2001 From: Javier Vegas-Regidor Date: Wed, 27 Oct 2021 07:02:51 +0200 Subject: [PATCH 11/11] progressing... --- esmvalcore/_recipe.py | 3 ++ .../cmor/_fixes/cordex/test_generic_cordex.py | 46 +++++++++++++++++++ 2 files changed, 49 insertions(+) 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/tests/integration/cmor/_fixes/cordex/test_generic_cordex.py b/tests/integration/cmor/_fixes/cordex/test_generic_cordex.py index ef3954518b..63a35dd177 100644 --- a/tests/integration/cmor/_fixes/cordex/test_generic_cordex.py +++ b/tests/integration/cmor/_fixes/cordex/test_generic_cordex.py @@ -1,4 +1,7 @@ """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 @@ -7,3 +10,46 @@ 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()