diff --git a/esmvalcore/cmor/tables/custom/CMOR_soz.dat b/esmvalcore/cmor/tables/custom/CMOR_soz.dat new file mode 100644 index 0000000000..725e454e16 --- /dev/null +++ b/esmvalcore/cmor/tables/custom/CMOR_soz.dat @@ -0,0 +1,22 @@ +!============ +variable_entry: soz +!============ +modeling_realm: atmos +!---------------------------------- +! Variable attributes: +!---------------------------------- +standard_name: equivalent_thickness_at_stp_of_atmosphere_ozone_content +units: m +cell_methods: time: mean +cell_measures: area: areacella +long_name: Stratospheric Ozone Column (O3 mole fraction >= 125 ppb) +comment: stratospheric ozone column calculated at 0 degrees C and 1 bar, such that 1m = 1e5 DU. Here, the stratosphere is defined as the region where O3 mole fraction >= 125 ppb. +!---------------------------------- +! Additional variable information: +!---------------------------------- +dimensions: longitude latitude time +type: real +valid_min: 0.0 +valid_max: 5000.0 +!---------------------------------- +! diff --git a/esmvalcore/cmor/tables/custom/CMOR_toz.dat b/esmvalcore/cmor/tables/custom/CMOR_toz.dat index b875dcbe57..d2de911497 100644 --- a/esmvalcore/cmor/tables/custom/CMOR_toz.dat +++ b/esmvalcore/cmor/tables/custom/CMOR_toz.dat @@ -1,4 +1,4 @@ -SOURCE: CCMI1 +SOURCE: CMIP6 !============ variable_entry: toz !============ @@ -6,12 +6,12 @@ modeling_realm: atmos !---------------------------------- ! Variable attributes: !---------------------------------- -standard_name: -units: DU +standard_name: equivalent_thickness_at_stp_of_atmosphere_ozone_content +units: m cell_methods: time: mean cell_measures: area: areacella long_name: Total Ozone Column -comment: total ozone column in DU +comment: Total ozone column calculated at 0 degrees C and 1 bar, such that 1m = 1e5 DU. !---------------------------------- ! Additional variable information: !---------------------------------- diff --git a/esmvalcore/cmor/tables/custom/CMOR_tozStderr.dat b/esmvalcore/cmor/tables/custom/CMOR_tozStderr.dat index 7d8769bc7c..0b80052bee 100644 --- a/esmvalcore/cmor/tables/custom/CMOR_tozStderr.dat +++ b/esmvalcore/cmor/tables/custom/CMOR_tozStderr.dat @@ -1,4 +1,4 @@ -SOURCE: CCMI1 +SOURCE: CMIP6 !============ variable_entry: tozStderr !============ @@ -6,12 +6,12 @@ modeling_realm: atmos !---------------------------------- ! Variable attributes: !---------------------------------- -standard_name: -units: DU +standard_name: equivalent_thickness_at_stp_of_atmosphere_ozone_content +units: m cell_methods: time: mean cell_measures: area: areacella long_name: Total Ozone Column Error -comment: total ozone column in DU +comment: Total ozone column error calculated at 0 degrees C and 1 bar, such that 1m = 1e5 DU. !---------------------------------- ! Additional variable information: !---------------------------------- diff --git a/esmvalcore/cmor/tables/custom/CMOR_troz.dat b/esmvalcore/cmor/tables/custom/CMOR_troz.dat new file mode 100644 index 0000000000..ea00615131 --- /dev/null +++ b/esmvalcore/cmor/tables/custom/CMOR_troz.dat @@ -0,0 +1,22 @@ +!============ +variable_entry: troz +!============ +modeling_realm: atmos +!---------------------------------- +! Variable attributes: +!---------------------------------- +standard_name: equivalent_thickness_at_stp_of_atmosphere_ozone_content +units: m +cell_methods: time: mean +cell_measures: area: areacella +long_name: Tropospheric Ozone Column (O3 mole fraction < 125 ppb) +comment: tropospheric ozone column calculated at 0 degrees C and 1 bar, such that 1m = 1e5 DU. Here, the troposphere is defined as the region where O3 mole fraction < 125 ppb. +!---------------------------------- +! Additional variable information: +!---------------------------------- +dimensions: longitude latitude time +type: real +valid_min: 0.0 +valid_max: 5000.0 +!---------------------------------- +! diff --git a/esmvalcore/preprocessor/_derive/_shared.py b/esmvalcore/preprocessor/_derive/_shared.py index e6d07011f1..d7f74274bc 100644 --- a/esmvalcore/preprocessor/_derive/_shared.py +++ b/esmvalcore/preprocessor/_derive/_shared.py @@ -170,15 +170,15 @@ def _create_pressure_array(cube, ps_cube, top_limit): def _get_pressure_level_widths(array, air_pressure_axis=1): """Compute pressure level widths. - For a 1D array with pressure level columns, return a 1D array with - pressure level widths. + For array with pressure level columns, return array with pressure + level widths. """ array = np.copy(array) if np.any(np.diff(array, axis=air_pressure_axis) > 0.0): raise ValueError("Pressure level value increased with height") - # Calculate centers + # Calculate array of centers between two neighboring pressure levels indices = [slice(None)] * array.ndim array_shifted = np.roll(array, -1, axis=air_pressure_axis) index_0 = deepcopy(indices) diff --git a/esmvalcore/preprocessor/_derive/soz.py b/esmvalcore/preprocessor/_derive/soz.py new file mode 100644 index 0000000000..1454ec505e --- /dev/null +++ b/esmvalcore/preprocessor/_derive/soz.py @@ -0,0 +1,107 @@ +"""Derivation of variable ``soz``.""" + +import dask.array as da +import iris + +from ._baseclass import DerivedVariableBase +from .toz import DerivedVariable as Toz +from .toz import add_longitude_coord, interpolate_hybrid_plevs + +# O3 mole fraction threshold (in ppb) that is used for the definition of the +# stratosphere (stratosphere = region where O3 mole fraction is at least as +# high as the threshold value) +STRATOSPHERIC_O3_THRESHOLD = 125.0 + + +class DerivedVariable(DerivedVariableBase): + """Derivation of variable ``soz``.""" + + @staticmethod + def required(project): + """Declare the variables needed for derivation.""" + if project == 'CMIP6': + required = [{'short_name': 'o3'}] + else: + required = [{'short_name': 'tro3'}] + return required + + @staticmethod + def calculate(cubes): + """Compute stratospheric column ozone. + + Note + ---- + Here, the stratosphere is defined as the region in which the O3 mole + fraction is at least as high as the given threshold + (``STRATOSPHERIC_O3_THRESHOLD``). + + In the calculation of ``toz``, the surface air pressure (``ps``) is + used to determine the pressure level width of the lowest layer. For + ``soz``, this lowest layer can be ignored since it is not located in + the stratosphere (it will be masked out due to the O3 mole fraction + threshold). Thus, the surface air pressure (``ps``) is not necessary + for the derivation of ``soz`` and is simply replaced with the lowest + pressure level in the data to be able to use the ``toz`` derivation + function. + + The calculation of ``soz`` consists of three steps: + (1) Mask out O3 mole fractions smaller than given threshold. + (2) Cut out the lowest pressure level from the data and use it as + surface air pressure (``toz``). + (3) Use derivation function of ``toz`` to calculate ``soz`` (using the + masked data). + + """ + o3_cube = cubes.extract_cube( + iris.Constraint(name='mole_fraction_of_ozone_in_air') + ) + + # If o3 is given on hybrid pressure levels (e.g., from Table AERmon), + # interpolate it to regular pressure levels + if len(o3_cube.coord_dims('air_pressure')) > 1: + o3_cube = interpolate_hybrid_plevs(o3_cube) + + # To support zonal mean o3 (e.g., from Table AERmonZ), add longitude + # coordinate if necessary + if not o3_cube.coords('longitude'): + o3_cube = add_longitude_coord(o3_cube) + + # (1) Mask O3 mole fraction using the given threshold + o3_cube.convert_units('1e-9') + mask = o3_cube.lazy_data() < STRATOSPHERIC_O3_THRESHOLD + mask |= da.ma.getmaskarray(o3_cube.lazy_data()) + o3_cube.data = da.ma.masked_array(o3_cube.lazy_data(), mask=mask) + + # (2) Add surrogate for the surface air pressure (ps) cube using the + # lowest pressure level available in the data (this is fine since the + # the lowest pressure level is far away from the stratosphere). + + # Get dummy ps cube with correct dimensions + ps_dims = (o3_cube.coord_dims('time') + + o3_cube.coord_dims('latitude') + + o3_cube.coord_dims('longitude')) + idx_to_extract_ps = [0] * o3_cube.ndim + for ps_dim in ps_dims: + idx_to_extract_ps[ps_dim] = slice(None) + ps_cube = o3_cube[tuple(idx_to_extract_ps)].copy() + + # Set ps data using lowest pressure level available and add correct + # metadata + lowest_plev = o3_cube.coord('air_pressure').points.max() + ps_data = da.broadcast_to(lowest_plev, ps_cube.shape) + ps_cube.data = ps_data + ps_cube.var_name = 'ps' + ps_cube.standard_name = 'surface_air_pressure' + ps_cube.long_name = 'Surface Air Pressure' + ps_cube.units = o3_cube.coord('air_pressure').units + + # Cut lowest pressure level from o3_cube + z_dim = o3_cube.coord_dims('air_pressure')[0] + idx_to_cut_lowest_plev = [slice(None)] * o3_cube.ndim + idx_to_cut_lowest_plev[z_dim] = slice(1, None) + o3_cube = o3_cube[tuple(idx_to_cut_lowest_plev)] + + # (3) Use derivation function of toz to calculate soz using the masked + # o3 cube and the surrogate ps cube + cubes = iris.cube.CubeList([o3_cube, ps_cube]) + return Toz.calculate(cubes) diff --git a/esmvalcore/preprocessor/_derive/toz.py b/esmvalcore/preprocessor/_derive/toz.py index 32fd9e8334..8590702ce4 100644 --- a/esmvalcore/preprocessor/_derive/toz.py +++ b/esmvalcore/preprocessor/_derive/toz.py @@ -1,9 +1,14 @@ -"""Derivation of variable `toz`.""" +"""Derivation of variable ``toz``.""" + +import warnings import cf_units import iris from scipy import constants +from esmvalcore.cmor.table import CMOR_TABLES + +from .._regrid import extract_levels, regrid from ._baseclass import DerivedVariableBase from ._shared import pressure_level_widths @@ -19,16 +24,53 @@ DOBSON_UNIT = cf_units.Unit('2.69e20 m^-2') +def add_longitude_coord(cube, ps_cube=None): + """Add dimensional ``longitude`` coordinate of length 1 to cube.""" + lon_coord = iris.coords.DimCoord([180.0], bounds=[[0.0, 360.0]], + var_name='lon', + standard_name='longitude', + long_name='longitude', + units='degrees_east') + new_dim_coords = [(c, cube.coord_dims(c)) for c in cube.dim_coords] + new_dim_coords.append((lon_coord, cube.ndim)) + new_aux_coords = [(c, cube.coord_dims(c)) for c in cube.aux_coords] + new_cube = iris.cube.Cube( + cube.core_data()[..., None], + dim_coords_and_dims=new_dim_coords, + aux_coords_and_dims=new_aux_coords, + ) + new_cube.metadata = cube.metadata + return new_cube + + +def interpolate_hybrid_plevs(cube): + """Interpolate hybrid pressure levels.""" + # Use CMIP6's plev19 target levels (in Pa) + target_levels = CMOR_TABLES['CMIP6'].coords['plev19'].requested + cube.coord('air_pressure').convert_units('Pa') + cube = extract_levels(cube, target_levels, 'linear', + coordinate='air_pressure') + return cube + + class DerivedVariable(DerivedVariableBase): - """Derivation of variable `toz`.""" + """Derivation of variable ``toz``.""" @staticmethod def required(project): """Declare the variables needed for derivation.""" + # TODO: make get_required _derive/__init__.py use variables as argument + # and make this dependent on mip if project == 'CMIP6': - required = [{'short_name': 'o3'}, {'short_name': 'ps'}] + required = [ + {'short_name': 'o3'}, + {'short_name': 'ps', 'mip': 'Amon'}, + ] else: - required = [{'short_name': 'tro3'}, {'short_name': 'ps'}] + required = [ + {'short_name': 'tro3'}, + {'short_name': 'ps'}, + ] return required @staticmethod @@ -41,24 +83,53 @@ def calculate(cubes): upper integration bound of 0 Pa is used. """ - tro3_cube = cubes.extract_cube( - iris.Constraint(name='mole_fraction_of_ozone_in_air')) + o3_cube = cubes.extract_cube( + iris.Constraint(name='mole_fraction_of_ozone_in_air') + ) ps_cube = cubes.extract_cube( - iris.Constraint(name='surface_air_pressure')) + iris.Constraint(name='surface_air_pressure') + ) + + # If o3 is given on hybrid pressure levels (e.g., from Table AERmon), + # interpolate it to regular pressure levels + if len(o3_cube.coord_dims('air_pressure')) > 1: + o3_cube = interpolate_hybrid_plevs(o3_cube) + + # To support zonal mean o3 (e.g., from Table AERmonZ), add longitude + # coordinate and collapsed ps cube if necessary to ensure that they + # have correct shapes + if not o3_cube.coords('longitude'): + o3_cube = add_longitude_coord(o3_cube) + ps_cube = ps_cube.collapsed('longitude', iris.analysis.MEAN) + ps_cube.remove_coord('longitude') + ps_cube = add_longitude_coord(ps_cube) + + # If the horizontal dimensions of ps and o3 differ, regrid ps + # Note: regrid() checks if the regridding is really necessary before + # running the actual interpolation + ps_cube = regrid(ps_cube, o3_cube, 'linear') - p_layer_widths = pressure_level_widths(tro3_cube, + # Actual derivation of toz using o3 mole fraction and pressure level + # widths + p_layer_widths = pressure_level_widths(o3_cube, ps_cube, top_limit=0.0) - toz_cube = (tro3_cube * p_layer_widths / STANDARD_GRAVITY * MW_O3 / + toz_cube = (o3_cube * p_layer_widths / STANDARD_GRAVITY * MW_O3 / MW_AIR) - toz_cube = toz_cube.collapsed('air_pressure', iris.analysis.SUM) - toz_cube.units = (tro3_cube.units * p_layer_widths.units / + with warnings.catch_warnings(): + warnings.filterwarnings( + 'ignore', category=UserWarning, + message='Collapsing a non-contiguous coordinate') + toz_cube = toz_cube.collapsed('air_pressure', iris.analysis.SUM) + toz_cube.units = (o3_cube.units * p_layer_widths.units / STANDARD_GRAVITY_UNIT * MW_O3_UNIT / MW_AIR_UNIT) - # Convert from kg m^-2 to Dobson unit (2.69e20 m^-2 ) + # Convert from kg m^-2 to Dobson units DU (2.69e20 m^-2 ) and from + # DU to m (1 mm = 100 DU) toz_cube = toz_cube / MW_O3 * AVOGADRO_CONST toz_cube.units = toz_cube.units / MW_O3_UNIT * AVOGADRO_CONST_UNIT toz_cube.convert_units(DOBSON_UNIT) - toz_cube.units = 'DU' + toz_cube.data = toz_cube.core_data() * 1e-5 + toz_cube.units = 'm' return toz_cube diff --git a/esmvalcore/preprocessor/_derive/troz.py b/esmvalcore/preprocessor/_derive/troz.py new file mode 100644 index 0000000000..d9b23fa713 --- /dev/null +++ b/esmvalcore/preprocessor/_derive/troz.py @@ -0,0 +1,61 @@ +"""Derivation of variable ``troz``.""" + +import dask.array as da +import iris + +from ._baseclass import DerivedVariableBase +from .soz import STRATOSPHERIC_O3_THRESHOLD +from .toz import DerivedVariable as Toz +from .toz import add_longitude_coord, interpolate_hybrid_plevs + + +class DerivedVariable(DerivedVariableBase): + """Derivation of variable ``troz``.""" + + @staticmethod + def required(project): + """Declare the variables needed for derivation.""" + return Toz.required(project) + + @staticmethod + def calculate(cubes): + """Compute tropospheric column ozone. + + Note + ---- + Here, the troposphere is defined as the region in which the O3 mole + fraction is smaller than the given threshold + (``STRATOSPHERIC_O3_THRESHOLD``). + + """ + o3_cube = cubes.extract_cube( + iris.Constraint(name='mole_fraction_of_ozone_in_air') + ) + ps_cube = cubes.extract_cube( + iris.Constraint(name='surface_air_pressure') + ) + + # If o3 is given on hybrid pressure levels (e.g., from Table AERmon), + # interpolate it to regular pressure levels + if len(o3_cube.coord_dims('air_pressure')) > 1: + o3_cube = interpolate_hybrid_plevs(o3_cube) + + # To support zonal mean o3 (e.g., from Table AERmonZ), add longitude + # coordinate and collapsed ps cube if necessary to ensure that they + # have correct shapes + if not o3_cube.coords('longitude'): + o3_cube = add_longitude_coord(o3_cube) + ps_cube = ps_cube.collapsed('longitude', iris.analysis.MEAN) + ps_cube.remove_coord('longitude') + ps_cube = add_longitude_coord(ps_cube) + + # Mask O3 mole fraction using the given threshold + o3_cube.convert_units('1e-9') + mask = o3_cube.lazy_data() >= STRATOSPHERIC_O3_THRESHOLD + mask |= da.ma.getmaskarray(o3_cube.lazy_data()) + o3_cube.data = da.ma.masked_array(o3_cube.lazy_data(), mask=mask) + + # Use derivation function of toz to calculate troz using the masked o3 + # cube and the ps cube + cubes = iris.cube.CubeList([o3_cube, ps_cube]) + return Toz.calculate(cubes) diff --git a/tests/integration/cmor/test_table.py b/tests/integration/cmor/test_table.py index 76fc3374a9..0be2c60dc9 100644 --- a/tests/integration/cmor/test_table.py +++ b/tests/integration/cmor/test_table.py @@ -137,7 +137,7 @@ def test_omon_ta_succes_if_strict(self): self.assertEqual(var.frequency, 'mon') def test_omon_toz_succes_if_strict(self): - """Get troz does not fail with Omon if not strict.""" + """Get toz does not fail with Omon if not strict.""" self.variables_info.strict = False var = self.variables_info.get_variable('Omon', 'toz') self.assertEqual(var.short_name, 'toz') @@ -296,7 +296,7 @@ def test_aermon_ta_succes_if_strict(self): self.assertEqual(var.frequency, 'mon') def test_omon_toz_succes_if_strict(self): - """Get troz does not fail with Omon if not strict.""" + """Get toz does not fail with Omon if not strict.""" self.variables_info.strict = False var = self.variables_info.get_variable('Omon', 'toz') self.assertEqual(var.short_name, 'toz') @@ -365,7 +365,7 @@ def test_aermon_ta_succes_if_strict(self): self.assertEqual(var.frequency, '') def test_omon_toz_succes_if_strict(self): - """Get troz does not fail with Omon if not strict.""" + """Get toz does not fail with Omon if not strict.""" self.variables_info.strict = False var = self.variables_info.get_variable('O1', 'toz') self.assertEqual(var.short_name, 'toz') diff --git a/tests/unit/preprocessor/_derive/test_co2s.py b/tests/unit/preprocessor/_derive/test_co2s.py index 3fd364edd7..1d6ff8ae4d 100644 --- a/tests/unit/preprocessor/_derive/test_co2s.py +++ b/tests/unit/preprocessor/_derive/test_co2s.py @@ -11,12 +11,14 @@ def get_coord_spec(include_plev=True): """Coordinate specs for cubes.""" time_coord = iris.coords.DimCoord([0], var_name='time', standard_name='time', - units='days since 0000-01-01 00:00:00') + units='days since 0001-01-01 00:00:00') lat_coord = iris.coords.DimCoord([0.0, 1.0], var_name='latitude', standard_name='latitude', units='degrees') lon_coord = iris.coords.DimCoord([0.0, 1.0], var_name='longitude', standard_name='longitude', units='degrees') + lat_coord.guess_bounds() + lon_coord.guess_bounds() if include_plev: plev_coord = iris.coords.DimCoord([100000.0, 90000.0, 50000.0], var_name='plev', diff --git a/tests/unit/preprocessor/_derive/test_soz.py b/tests/unit/preprocessor/_derive/test_soz.py new file mode 100644 index 0000000000..a523c2570e --- /dev/null +++ b/tests/unit/preprocessor/_derive/test_soz.py @@ -0,0 +1,117 @@ +"""Test derivation of ``soz``.""" +import dask.array as da +import iris +import numpy as np +import pytest + +import esmvalcore.preprocessor._derive.soz as soz + +from .test_toz import get_masked_o3_cube, get_masked_o3_hybrid_plevs_cube + + +def get_o3_cube(): + """Get ``o3`` input cube.""" + o3_cube = get_masked_o3_cube() + o3_cube.data = da.ma.masked_greater( + [[ + [[500.0, 700.0], [800.0, 900.0]], + [[1251.0, 1249.0], [1260.0, 1200.0]], + [[1000.0, 2000.0], [3000.0, 12000.0]], + ]], + 10000.0, + ) + o3_cube.units = '1e-10' + return o3_cube + + +@pytest.fixture +def cubes(): + """Input cubes for derivation of ``soz``.""" + o3_cube = get_o3_cube() + return iris.cube.CubeList([o3_cube]) + + +@pytest.fixture +def cubes_no_lon(): + """Zonal mean input cubes for derivation of ``soz``.""" + o3_cube = get_o3_cube() + o3_cube = o3_cube.collapsed('longitude', iris.analysis.MEAN) + o3_cube.remove_coord('longitude') + return iris.cube.CubeList([o3_cube]) + + +@pytest.fixture +def cubes_hybrid_plevs(): + """Input cubes with hybrid pressure levels for derivation of ``soz``.""" + o3_cube = get_masked_o3_hybrid_plevs_cube() + o3_cube.data = da.ma.masked_greater( + [[ + [[500.0, 700.0], [800.0, 900.0]], + [[1251.0, 1249.0], [1260.0, 1200.0]], + [[1000.0, 2000.0], [3000.0, 12000.0]], + ]], + 10000.0, + ) + o3_cube.units = '1e-10' + return iris.cube.CubeList([o3_cube]) + + +def test_soz_calculate(cubes): + """Test function ``calculate``.""" + derived_var = soz.DerivedVariable() + + out_cube = derived_var.calculate(cubes) + + assert out_cube.units == 'm' + assert out_cube.shape == (1, 2, 2) + expected_data = np.ma.masked_invalid( + [[ + [29.543266581831194e-5, 110.2066965482645e-5], + [195.06585289042815e-5, np.nan], + ]], + ) + expected_mask = [[[False, False], [False, True]]] + np.testing.assert_allclose(out_cube.data, expected_data) + np.testing.assert_allclose(out_cube.data.mask, expected_mask) + + +def test_soz_calculate_no_lon(cubes_no_lon): + """Test function ``calculate`` for zonal mean cubes.""" + derived_var = soz.DerivedVariable() + + out_cube = derived_var.calculate(cubes_no_lon) + + assert out_cube.units == 'm' + assert out_cube.shape == (1, 2, 1) + assert not np.ma.is_masked(out_cube.data) + np.testing.assert_allclose( + out_cube.data, [[[82.65502241119836e-5], [165.31004482239672e-5]]] + ) + + +def test_soz_calculate_hybrid_plevs(cubes_hybrid_plevs): + """Test function ``calculate`` for cubes with hybrid pressure levels.""" + derived_var = soz.DerivedVariable() + + out_cube = derived_var.calculate(cubes_hybrid_plevs) + + assert out_cube.units == 'm' + assert out_cube.shape == (1, 2, 2) + expected_data = np.ma.masked_invalid( + [[[np.nan, 32.40347475318536e-5], [44.53039332403313e-5, np.nan]]] + ) + expected_mask = [[[True, False], [False, True]]] + np.testing.assert_allclose(out_cube.data, expected_data) + np.testing.assert_allclose(out_cube.data.mask, expected_mask) + + +@pytest.mark.parametrize('project,out', [ + ('CMIP5', [{'short_name': 'tro3'}]), + ('TEST', [{'short_name': 'tro3'}]), + ('CMIP6', [{'short_name': 'o3'}]), +]) +def test_soz_required(project, out): + """Test function ``required``.""" + derived_var = soz.DerivedVariable() + output = derived_var.required(project) + assert output == out diff --git a/tests/unit/preprocessor/_derive/test_toz.py b/tests/unit/preprocessor/_derive/test_toz.py index 262deab1e0..ab794de4be 100644 --- a/tests/unit/preprocessor/_derive/test_toz.py +++ b/tests/unit/preprocessor/_derive/test_toz.py @@ -5,19 +5,21 @@ import pytest import esmvalcore.preprocessor._derive.toz as toz + from .test_co2s import get_coord_spec, get_ps_cube -@pytest.fixture -def masked_cubes(): - """Masked O3 cube.""" +def get_masked_o3_cube(): + """Get masked ``o3`` cube.""" coord_spec = get_coord_spec() - o3_data = da.ma.masked_less([[[[0.0, -1.0], - [-1.0, -1.0]], - [[1.0, 2.0], - [3.0, -1.0]], - [[2.0, 2.0], - [2.0, 2.0]]]], 0.0) + o3_data = da.ma.masked_less( + [[ + [[0.0, -1.0], [-1.0, -1.0]], + [[1.0, 2.0], [3.0, -1.0]], + [[2.0, 2.0], [2.0, 2.0]], + ]], + 0.0, + ) o3_cube = iris.cube.Cube( o3_data, var_name='o3', @@ -25,7 +27,59 @@ def masked_cubes(): units='1e-9', dim_coords_and_dims=coord_spec, ) + return o3_cube + + +def get_masked_o3_hybrid_plevs_cube(): + """Get masked ``o3`` cube with hybrid pressure levels.""" + o3_cube = get_masked_o3_cube() + o3_cube.remove_coord('air_pressure') + + ap_coord = iris.coords.AuxCoord([0.0, 10000.0, 0.0], var_name='ap', + units='Pa') + b_coord = iris.coords.AuxCoord([0.95, 0.8, 0.7], var_name='b', units='1') + ps_coord = iris.coords.AuxCoord( + [[[100000.0, 100000.0], [100000.0, 100000.0]]], var_name='ps', + units='Pa') + z_coord = iris.coords.DimCoord([0.95, 0.9, 0.7], var_name='lev', + units='1', attributes={'positive': 'down'}) + o3_cube.add_aux_coord(ap_coord, 1) + o3_cube.add_aux_coord(b_coord, 1) + o3_cube.add_aux_coord(ps_coord, (0, 2, 3)) + o3_cube.add_dim_coord(z_coord, 1) + + aux_factory = iris.aux_factory.HybridPressureFactory( + delta=ap_coord, sigma=b_coord, surface_air_pressure=ps_coord) + o3_cube.add_aux_factory(aux_factory) + + return o3_cube + + +@pytest.fixture +def masked_cubes(): + """Masked O3 cube.""" + o3_cube = get_masked_o3_cube() + ps_cube = get_ps_cube() + return iris.cube.CubeList([o3_cube, ps_cube]) + + +@pytest.fixture +def masked_cubes_no_lon(): + """Masked zonal mean O3 cube.""" + o3_cube = get_masked_o3_cube() + o3_cube = o3_cube.collapsed('longitude', iris.analysis.MEAN) + o3_cube.remove_coord('longitude') + ps_cube = get_ps_cube() + ps_cube.data = [[[101300.0, 101300.0], [101300.0, 101300.0]]] + return iris.cube.CubeList([o3_cube, ps_cube]) + + +@pytest.fixture +def masked_cubes_hybrid_plevs(): + """Masked zonal mean O3 cube on hybrid levels.""" + o3_cube = get_masked_o3_hybrid_plevs_cube() ps_cube = get_ps_cube() + ps_cube.data = [[[101300.0, 101300.0], [101300.0, 101300.0]]] return iris.cube.CubeList([o3_cube, ps_cube]) @@ -33,12 +87,13 @@ def masked_cubes(): def unmasked_cubes(): """Unmasked O3 cube.""" coord_spec = get_coord_spec() - o3_data = da.array([[[[2.0, 1.0], - [0.8, 1.0]], - [[1.5, 0.8], - [2.0, 3.0]], - [[4.0, 1.0], - [3.0, 2.0]]]]) + o3_data = da.array( + [[ + [[2.0, 1.0], [0.8, 1.0]], + [[1.5, 0.8], [2.0, 3.0]], + [[4.0, 1.0], [3.0, 2.0]], + ]], + ) o3_cube = iris.cube.Cube( o3_data, var_name='o3', @@ -53,29 +108,68 @@ def unmasked_cubes(): def test_toz_calculate_masked_cubes(masked_cubes): """Test function ``calculate`` with masked cube.""" derived_var = toz.DerivedVariable() + out_cube = derived_var.calculate(masked_cubes) + + assert out_cube.units == 'm' assert not np.ma.is_masked(out_cube.data) - np.testing.assert_allclose(out_cube.data, - [[[1.2988646378902597, 0.7871906896304607], - [1.6924599827054907, 0.9446288275565529]]]) - assert out_cube.units == 'DU' + np.testing.assert_allclose( + out_cube.data, + [[ + [1.2988646378902597e-5, 0.7871906896304607e-5], + [1.6924599827054907e-5, 0.9446288275565529e-5], + ]], + ) + + +def test_toz_calculate_masked_cubes_no_lon(masked_cubes_no_lon): + """Test function ``calculate`` with zonal mean masked cube.""" + derived_var = toz.DerivedVariable() + + out_cube = derived_var.calculate(masked_cubes_no_lon) + + assert out_cube.units == 'm' + assert not np.ma.is_masked(out_cube.data) + np.testing.assert_allclose( + out_cube.data, [[[1.3972634740940675e-5], [1.6924599827054907e-5]]], + ) + + +def test_toz_calculate_masked_cubes_hybrid_plevs(masked_cubes_hybrid_plevs): + """Test function ``calculate`` with zonal mean masked cube.""" + derived_var = toz.DerivedVariable() + + out_cube = derived_var.calculate(masked_cubes_hybrid_plevs) + + assert out_cube.units == 'm' + assert not np.ma.is_masked(out_cube.data) + np.testing.assert_allclose( + out_cube.data, + [[ + [0.33701601399804104e-5, 0.3739155775744688e-5], + [0.440334792012039e-5, 0.19679767240761517e-5], + ]], + ) def test_toz_calculate_unmasked_cubes(unmasked_cubes): """Test function ``calculate`` with unmasked cube.""" derived_var = toz.DerivedVariable() + out_cube = derived_var.calculate(unmasked_cubes) + + assert out_cube.units == 'm' assert not np.ma.is_masked(out_cube.data) - np.testing.assert_allclose(out_cube.data, - [[[2.65676858, 0.39359534], - [2.04669579, 0.94462883]]]) - assert out_cube.units == 'DU' + np.testing.assert_allclose( + out_cube.data, + [[[2.65676858e-5, 0.39359534e-5], [2.04669579e-5, 0.94462883e-5]]], + ) @pytest.mark.parametrize('project,out', [ ('CMIP5', [{'short_name': 'tro3'}, {'short_name': 'ps'}]), ('TEST', [{'short_name': 'tro3'}, {'short_name': 'ps'}]), - ('CMIP6', [{'short_name': 'o3'}, {'short_name': 'ps'}]), + ('CMIP6', [{'short_name': 'o3'}, {'short_name': 'ps', 'mip': 'Amon'}]), ]) def test_toz_required(project, out): """Test function ``required``.""" diff --git a/tests/unit/preprocessor/_derive/test_troz.py b/tests/unit/preprocessor/_derive/test_troz.py new file mode 100644 index 0000000000..e31bf83b1f --- /dev/null +++ b/tests/unit/preprocessor/_derive/test_troz.py @@ -0,0 +1,127 @@ +"""Test derivation of ``troz``.""" +import iris +import numpy as np +import pytest +from iris.coords import AuxCoord, DimCoord +from iris.util import broadcast_to_shape + +import esmvalcore.preprocessor._derive.troz as troz + +from .test_toz import get_masked_o3_cube, get_ps_cube + + +def get_o3_cube(): + """Get ``o3`` input cube.""" + o3_cube = get_masked_o3_cube() + o3_cube.data = [[ + [[50.0, 70.0], [80.0, 90.0]], + [[70.0, 90.0], [100.0, 110.0]], + [[130, 140.0], [150.0, 160.0]], + ]] + o3_cube.units = '1e-9' + return o3_cube + + +@pytest.fixture +def cubes(): + """Input cubes for derivation of ``troz``.""" + o3_cube = get_o3_cube() + ps_cube = get_ps_cube() + ps_cube.data = [[[101300.0, 101300.0], [101300.0, 101300.0]]] + return iris.cube.CubeList([o3_cube, ps_cube]) + + +@pytest.fixture +def cubes_no_lon(): + """Zonal mean input cubes for derivation of ``troz``.""" + o3_cube = get_o3_cube() + o3_cube = o3_cube.collapsed('longitude', iris.analysis.MEAN) + o3_cube.remove_coord('longitude') + ps_cube = get_ps_cube() + ps_cube.data = [[[101300.0, 101300.0], [101300.0, 101300.0]]] + return iris.cube.CubeList([o3_cube, ps_cube]) + + +@pytest.fixture +def cubes_hybrid_plevs(): + """Input cubes with hybrid pressure levels for derivation of ``troz``.""" + o3_cube = get_o3_cube() + plev_coord = o3_cube.coord('air_pressure') + hybrid_plev_coord = AuxCoord( + broadcast_to_shape( + plev_coord.points, o3_cube.shape, o3_cube.coord_dims(plev_coord) + ), + ) + hybrid_plev_coord.metadata = plev_coord.metadata + alt_coord = DimCoord( + [0.0, 1000.0, 3000.0], + standard_name='altitude', + attributes={'positive': 'up'}, + ) + o3_cube.remove_coord(plev_coord) + o3_cube.add_aux_coord(hybrid_plev_coord, (0, 1, 2, 3)) + o3_cube.add_dim_coord(alt_coord, 1) + + ps_cube = get_ps_cube() + ps_cube.data = [[[101300.0, 101300.0], [101300.0, 101300.0]]] + + return iris.cube.CubeList([o3_cube, ps_cube]) + + +def test_troz_calculate(cubes): + """Test function ``calculate``.""" + derived_var = troz.DerivedVariable() + + out_cube = derived_var.calculate(cubes) + + assert out_cube.units == 'm' + assert out_cube.shape == (1, 2, 2) + assert not np.ma.is_masked(out_cube.data) + expected_data = [[ + [16.255487740869038e-5, 21.1833014579557e-5], + [23.647208316499057e-5, 26.111115175042404e-5], + ]] + np.testing.assert_allclose(out_cube.data, expected_data) + + +def test_troz_calculate_no_lon(cubes_no_lon): + """Test function ``calculate`` for zonal mean cubes.""" + derived_var = troz.DerivedVariable() + + out_cube = derived_var.calculate(cubes_no_lon) + + assert out_cube.units == 'm' + assert out_cube.shape == (1, 2, 1) + assert not np.ma.is_masked(out_cube.data) + print(out_cube.data) + np.testing.assert_allclose( + out_cube.data, [[[18.71939459941235e-5], [24.87916174577070e-5]]] + ) + + +def test_troz_calculate_hybrid_plevs(cubes_hybrid_plevs): + """Test function ``calculate`` for cubes with hybrid pressure levels.""" + derived_var = troz.DerivedVariable() + + out_cube = derived_var.calculate(cubes_hybrid_plevs) + + assert out_cube.units == 'm' + assert out_cube.shape == (1, 2, 2) + assert not np.ma.is_masked(out_cube.data) + expected_data = [[ + [31.581106479612044e-5, 27.640725083830575e-5], + [18.198372761713192e-5, 20.071886603033692e-5], + ]] + np.testing.assert_allclose(out_cube.data, expected_data) + + +@pytest.mark.parametrize('project,out', [ + ('CMIP5', [{'short_name': 'tro3'}, {'short_name': 'ps'}]), + ('TEST', [{'short_name': 'tro3'}, {'short_name': 'ps'}]), + ('CMIP6', [{'short_name': 'o3'}, {'short_name': 'ps', 'mip': 'Amon'}]), +]) +def test_toz_required(project, out): + """Test function ``required``.""" + derived_var = troz.DerivedVariable() + output = derived_var.required(project) + assert output == out