From 09b25a6a681b7f9d47c89970f7bd096ad825d5df Mon Sep 17 00:00:00 2001 From: Javier Vegas-Regidor Date: Wed, 6 Oct 2021 08:25:21 +0200 Subject: [PATCH 1/2] Fix interpolation bug --- earthdiagnostics/ocean/interpolatecdo.py | 134 ++++++----------------- environment.yml | 2 +- test/unit/ocean/test_interpolatecdo.py | 61 +---------- 3 files changed, 35 insertions(+), 162 deletions(-) diff --git a/earthdiagnostics/ocean/interpolatecdo.py b/earthdiagnostics/ocean/interpolatecdo.py index c12632f2..36fa9b32 100644 --- a/earthdiagnostics/ocean/interpolatecdo.py +++ b/earthdiagnostics/ocean/interpolatecdo.py @@ -4,6 +4,7 @@ import os from cdo import CDOException import numpy as np +import iris from earthdiagnostics.constants import Basins from earthdiagnostics.diagnostic import ( @@ -141,7 +142,6 @@ class InterpolateCDO(Diagnostic): ), DiagnosticBoolOption("mask_oceans", True), DiagnosticOption("original_grid", ""), - DiagnosticBoolOption("weights_from_mask", True), DiagnosticFrequencyOption( default_value=diags.config.frequency), ) @@ -153,31 +153,24 @@ class InterpolateCDO(Diagnostic): raise Exception("Target grid not provided") job_list = list() weights = TempFile.get() - method = options["method"].lower() - if options["weights_from_mask"]: - temp = cls.get_sample_grid_file() - cls.compute_weights(method, target_grid, temp, weights) - os.remove(temp) - weights_job = None - else: - ( - startdate, - member, - chunk, - ) = diags.config.experiment.get_chunk_list()[0] - weights_job = ComputeWeights( - diags.data_manager, - startdate, - member, - chunk, - options["domain"], - options["variables"][0], - target_grid, - options["original_grid"], - weights, - options["method"], - options["frequency"], - ) + ( + startdate, + member, + chunk, + ) = diags.config.experiment.get_chunk_list()[0] + weights_job = ComputeWeights( + diags.data_manager, + startdate, + member, + chunk, + options["domain"], + options["variables"][0], + target_grid, + options["original_grid"], + weights, + options["method"], + options["frequency"], + ) for var in options["variables"]: for ( @@ -221,6 +214,7 @@ class InterpolateCDO(Diagnostic): Path to the file to store the weights """ + sample_file = cls.remove_redundant_column(sample_file) if method == InterpolateCDO.BILINEAR: Utils.cdo().genbil(target_grid, input=sample_file, output=weights) elif method == InterpolateCDO.BICUBIC: @@ -230,81 +224,6 @@ class InterpolateCDO(Diagnostic): elif method == InterpolateCDO.CONSERVATIVE2: Utils.cdo().gencon2(target_grid, input=sample_file, output=weights) - @classmethod - def get_sample_grid_file(cls): - """ - Get a sample grid file - - Create a sample grid file from the definition in the masks file - - Returns - ------- - str - """ - temp = TempFile.get() - - handler = Utils.open_cdf("mask.nc") - - lat_name, lon_name = cls._get_lat_lon_alias(handler) - lon_bnds_name = "{0}_bnds".format(lon_name) - lat_bnds_name = "{0}_bnds".format(lat_name) - - Utils.nco().ncks( - input="mask.nc", - output=temp, - options=( - "-O -v tmask,{0},{1},gphif,glamf".format(lat_name, lon_name), - ), - ) - handler = Utils.open_cdf(temp) - lon = handler.variables[lon_name] - lon.units = "degrees_east" - lon.long_name = "Longitude" - lon.nav_model = "Default grid" - lon.standard_name = "longitude" - lon.short_name = lon_name - lon.bounds = lon_bnds_name - - lat = handler.variables[lat_name] - lat.units = "degrees_north" - lat.long_name = "Latitude" - lat.nav_model = "Default grid" - lat.standard_name = "latitude" - lat.short_name = lat_name - lat.bounds = lat_bnds_name - - handler.createDimension("bounds", 4) - - lon_bnds = handler.createVariable( - lon_bnds_name, lon.datatype, ("j", "i", "bounds") - ) - corner_lat = handler.variables["glamf"][0, ...] - lon_bnds[:, :, 0] = corner_lat - lon_bnds[:, :, 1] = np.roll(corner_lat, 1, 0) - lon_bnds[:, :, 2] = np.roll(corner_lat, -1, 1) - lon_bnds[:, :, 3] = np.roll(lon_bnds[:, :, 1], -1, 1) - - lat_bnds = handler.createVariable( - lat_bnds_name, lat.datatype, ("j", "i", "bounds") - ) - corner_lat = handler.variables["gphif"][0, ...] - lat_bnds[:, :, 0] = corner_lat - lat_bnds[:, :, 1] = np.roll(corner_lat, 1, 0) - lat_bnds[:, :, 2] = np.roll(corner_lat, 1, 1) - lat_bnds[:, :, 3] = np.roll(lat_bnds[:, :, 1], 1, 1) - lat_bnds[0, :, 1] = lat_bnds[1, 0, 1] - 1 - lat_bnds[0, :, 3] = lat_bnds[1, 0, 3] - 1 - - tmask = handler.variables["tmask"] - tmask.coordinates = "time lev {0} {1}".format(lat_name, lon_name) - - handler.close() - - Utils.nco().ncks( - input=temp, output=temp, options=("-O -x -v gphif,glamf",) - ) - return temp - @classmethod def _get_lat_lon_alias(cls, handler): lat_name = None @@ -401,6 +320,8 @@ class InterpolateCDO(Diagnostic): var[:] = mask * var[:] handler.close() + variable_file = self.remove_redundant_column(variable_file) + temp = TempFile.get() try: Utils.cdo().remap( @@ -434,6 +355,17 @@ class InterpolateCDO(Diagnostic): self.regridded.set_local_file(temp) + @staticmethod + def remove_redundant_column(variable_file): + cube = iris.load_cube(variable_file) + temp = TempFile.get() + Utils.cdo().selindexbox( + 1, cube.shape[-1] - 1, 1, cube.shape[-2], + input=variable_file, + output=temp + ) + return temp + class ComputeWeights(Diagnostic): """ diff --git a/environment.yml b/environment.yml index 2808ec2d..521a6c73 100644 --- a/environment.yml +++ b/environment.yml @@ -6,7 +6,7 @@ channels: dependencies: - python<3.9 -- cdo +- cdo=1.9.8 - nco - eccodes - six diff --git a/test/unit/ocean/test_interpolatecdo.py b/test/unit/ocean/test_interpolatecdo.py index 5b546dd9..4170d32f 100644 --- a/test/unit/ocean/test_interpolatecdo.py +++ b/test/unit/ocean/test_interpolatecdo.py @@ -31,20 +31,12 @@ class TestInterpolate(TestCase): raise DiagnosticOptionError return value.split("-") - @patch( - "earthdiagnostics.ocean.interpolatecdo.InterpolateCDO.compute_weights" - ) - @patch( - "earthdiagnostics.ocean.interpolatecdo.InterpolateCDO." - "get_sample_grid_file" - ) @patch.object(DiagnosticVariableListOption, "parse", fake_parse) @patch("os.remove") @patch("earthdiagnostics.utils.TempFile.get") def test_generate_jobs( - self, mock_weights, mock_grid_file, mock_remove, mock_get + self, mock_remove, mock_get ): - mock_weights.return_value = None mock_get.return_value = "path_to_weights" jobs = InterpolateCDO.generate_jobs( @@ -262,56 +254,6 @@ class TestInterpolate(TestCase): "bicubic", "false", "orig", - "false", - ], - ) - self.assertEqual(len(jobs), 2) - self.assertEqual( - jobs[0], - InterpolateCDO( - self.data_manager, - "20010101", - 0, - 0, - ModelingRealms.ocean, - "var", - "target_grid", - "model_version", - False, - "orig", - None, - Frequencies.daily, - ), - ) - self.assertEqual( - jobs[1], - InterpolateCDO( - self.data_manager, - "20010101", - 0, - 1, - ModelingRealms.ocean, - "var", - "target_grid", - "model_version", - False, - "orig", - None, - Frequencies.daily, - ), - ) - - jobs = InterpolateCDO.generate_jobs( - self.diags, - [ - "interpcdo", - "ocean", - "var", - "target_grid", - "bicubic", - "false", - "orig", - "false", "mon" ], ) @@ -364,7 +306,6 @@ class TestInterpolate(TestCase): "bicubic", "false", "orig", - "false", "daily," "extra", ], -- GitLab From bf08c0a78fb272fb2b8bf0274610958e224ef00d Mon Sep 17 00:00:00 2001 From: Javier Vegas-Regidor Date: Wed, 6 Oct 2021 08:33:05 +0200 Subject: [PATCH 2/2] Add extra check --- earthdiagnostics/ocean/interpolatecdo.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/earthdiagnostics/ocean/interpolatecdo.py b/earthdiagnostics/ocean/interpolatecdo.py index 36fa9b32..9873f8c8 100644 --- a/earthdiagnostics/ocean/interpolatecdo.py +++ b/earthdiagnostics/ocean/interpolatecdo.py @@ -358,6 +358,8 @@ class InterpolateCDO(Diagnostic): @staticmethod def remove_redundant_column(variable_file): cube = iris.load_cube(variable_file) + if cube.coord('longitude').ndim == 1: + return variable_file temp = TempFile.get() Utils.cdo().selindexbox( 1, cube.shape[-1] - 1, 1, cube.shape[-2], -- GitLab