From bc2ae2ada5e29446fe0065ed554ed2da1ef40e45 Mon Sep 17 00:00:00 2001 From: aperez Date: Tue, 23 Jul 2019 11:25:37 +0200 Subject: [PATCH] Adapt regridding to CORDEX needs --- esmvalcore/preprocessor/_regrid.py | 15 +++++++++ esmvalcore/preprocessor/_regrid_esmpy.py | 40 ++++++++++++++---------- 2 files changed, 39 insertions(+), 16 deletions(-) diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py index 58a96491d4..4c88596c05 100644 --- a/esmvalcore/preprocessor/_regrid.py +++ b/esmvalcore/preprocessor/_regrid.py @@ -8,6 +8,7 @@ import numpy as np import stratify from iris.analysis import AreaWeighted, Linear, Nearest, UnstructuredNearest +from iris import coord_systems from ..cmor.fix import fix_file, fix_metadata from ..cmor.table import CMOR_TABLES @@ -250,6 +251,20 @@ def regrid(cube, target_grid, scheme, lat_offset=True, lon_offset=True): [coord] = coords cube.remove_coord(coord) + # we first define the name of the coordinate system for the target grid in order to avoid iris problems with the + # coordinate system given by WGS84. + if target_grid.coord_system() == None: + target_grid.coord('latitude').coord_system = iris.coord_systems.GeogCS(semi_major_axis=6378137.0, + semi_minor_axis=6356752.31424) + target_grid.coord('longitude').coord_system = iris.coord_systems.GeogCS(semi_major_axis=6378137.0, + semi_minor_axis=6356752.31424) + + if cube.coord_system() == None: + cube.coord('latitude').coord_system = iris.coord_systems.GeogCS(semi_major_axis=6378137.0, + semi_minor_axis=6356752.31424) + cube.coord('longitude').coord_system = iris.coord_systems.GeogCS(semi_major_axis=6378137.0, + semi_minor_axis=6356752.31424) + # Perform the horizontal regridding. if _attempt_irregular_regridding(cube, scheme): cube = esmpy_regrid(cube, target_grid, scheme) diff --git a/esmvalcore/preprocessor/_regrid_esmpy.py b/esmvalcore/preprocessor/_regrid_esmpy.py index 4ac6192451..3ec54bb140 100755 --- a/esmvalcore/preprocessor/_regrid_esmpy.py +++ b/esmvalcore/preprocessor/_regrid_esmpy.py @@ -6,17 +6,18 @@ import numpy as np from ._mapping import get_empty_data, map_slices, ref_to_dims_index +from esmvalcore.cmor._fixes.CORDEX.Cordex_boundaries_guessing import CordexFix ESMF_MANAGER = ESMF.Manager(debug=False) ESMF_LON, ESMF_LAT = 0, 1 -ESMF_REGRID_METHODS = { - 'linear': ESMF.RegridMethod.BILINEAR, - 'area_weighted': ESMF.RegridMethod.CONSERVE, - 'nearest': ESMF.RegridMethod.NEAREST_STOD, -} +#ESMF_REGRID_METHODS = { +# 'linear': ESMF.RegridMethod.BILINEAR, +# 'area_weighted': ESMF.RegridMethod.CONSERVE, +# 'nearest': ESMF.RegridMethod.NEAREST_STOD, +#} MASK_REGRIDDING_MASK_VALUE = { ESMF.RegridMethod.BILINEAR: np.array([1]), @@ -24,13 +25,13 @@ ESMF.RegridMethod.NEAREST_STOD: np.array([]), } -# ESMF_REGRID_METHODS = { -# 'bilinear': ESMF.RegridMethod.BILINEAR, -# 'patch': ESMF.RegridMethod.PATCH, -# 'conserve': ESMF.RegridMethod.CONSERVE, -# 'nearest_stod': ESMF.RegridMethod.NEAREST_STOD, -# 'nearest_dtos': ESMF.RegridMethod.NEAREST_DTOS, -# } +ESMF_REGRID_METHODS = { + 'bilinear': ESMF.RegridMethod.BILINEAR, + 'patch': ESMF.RegridMethod.PATCH, + 'conserve': ESMF.RegridMethod.CONSERVE, + 'nearest_stod': ESMF.RegridMethod.NEAREST_STOD, + 'nearest_dtos': ESMF.RegridMethod.NEAREST_DTOS, +} def cf_2d_bounds_to_esmpy_corners(bounds, circular): @@ -101,7 +102,7 @@ def get_grid(esmpy_lat, esmpy_lon, return grid -def is_lon_circular(lon): +def is_lon_circular(lon, cube): """Determine if longitudes are circular.""" if isinstance(lon, iris.coords.DimCoord): circular = lon.circular @@ -109,8 +110,15 @@ def is_lon_circular(lon): if lon.ndim == 1: seam = lon.bounds[-1, 1] - lon.bounds[0, 0] elif lon.ndim == 2: - seam = (lon.bounds[1:-1, -1, (1, 2)] - - lon.bounds[1:-1, 0, (0, 3)]) + try: + seam = (lon.bounds[1:-1, -1, (1, 2)] - lon.bounds[1:-1, 0, (0, 3)]) + + except: + object = CordexFix() + cube = object.fix_data(cube) + lon = cube.coord('longitude') + seam = (lon.bounds[1:-1, -1, (1, 2)] - lon.bounds[1:-1, 0, (0, 3)]) + else: raise NotImplementedError('AuxCoord longitude is higher ' 'dimensional than 2d. Giving up.') @@ -125,7 +133,7 @@ def cube_to_empty_field(cube): """Build an empty ESMF field from a cube.""" lat = cube.coord('latitude') lon = cube.coord('longitude') - circular = is_lon_circular(lon) + circular = is_lon_circular(lon, cube) esmpy_coords = coords_iris_to_esmpy(lat, lon, circular) grid = get_grid(*esmpy_coords, circular=circular) field = ESMF.Field(grid,