Added land/sea fraction weighting preprocessor

doc/esmvalcore/preprocessor.rst

@@ -10,6 +10,7 @@ following the default order in which they are applied:
 * :ref:`Variable derivation`
 * :ref:`CMOR check and dataset-specific fixes`
 * :ref:`Vertical interpolation`
+* :ref:`Weighting`
 * :ref:`Land/Sea/Ice masking`
 * :ref:`Horizontal regridding`
 * :ref:`Masking of missing values`
@@ -244,6 +245,50 @@ extract the levels and vertically regrid onto the vertical levels of
 	  points will be masked, otherwise they will be set to NaN.
 
 
+.. _weighting:
+
+Weighting
+=========
+
+.. _land/sea fraction weighting:
+
+Land/sea fraction weighting
+---------------------------
+
+This preprocessor allows weighting of data by land or sea fractions. In other
+words, this function multiplies the given input field by a fraction in the range 0-1 to
+account for the fact that not all grid points are completely land- or sea-covered.
+
+The application of this preprocessor is very important for most carbon cycle variables (and
+other land surface outputs), which are e.g. reported in units of
+:math:`kgC~m^{-2}`. Here, the surface unit actually refers to 'square meter of land/sea' and
+NOT 'square meter of gridbox'. In order to integrate these globally or
+regionally one has to weight by both the surface quantity and the
+land/sea fraction.
+
+For example, to weight an input field with the land fraction, the following
+preprocessor can be used:
+
+.. code-block:: yaml
+
+    preprocessors:
+      preproc_weighting:
+        weighting_landsea_fraction:
+          area_type: land
+          strict: true
+
+Allowed arguments for the keyword ``area_type`` are ``land`` (fraction is 1
+for grid cells with only land surface, 0 for grid cells with only sea surface
+and values in between 0 and 1 for coastal regions) and ``sea`` (1 for
+sea, 0 for land, in between for coastal regions). The optional argument
+``strict`` defines the behavior of this function in case the weighting was not
+possible (due to missing data or incompatible shapes). In the case of
+``strict: true`` (default), it fails, in the case of ``strict: false``, it does
+not apply the weighting and prints a debug message.
+
+See also :func:`esmvalcore.preprocessor.weighting_landsea_fraction`.
+
+
 .. _masking:
 
 Masking

esmvalcore/_recipe.py

@@ -5,6 +5,7 @@
 import re
 from collections import OrderedDict
 from copy import deepcopy
+from pprint import pformat
 
 import yaml
 from netCDF4 import Dataset
@@ -379,7 +380,6 @@ def _add_fxvar_keys(fx_var_dict, variable):
         fx_variable['mip'] = 'fx'
     # add missing cmor info
     _add_cmor_info(fx_variable, override=True)
-
     return fx_variable
 
 
@@ -408,28 +408,29 @@ def _get_correct_fx_file(variable, fx_varname, config_user):
     if fx_files:
         fx_files = fx_files[0]
 
-    logger.info("Using fx files for masking for %s of dataset %s:\n%s",
-                fx_var['short_name'], fx_var['dataset'],
-                fx_files)
-
     return fx_files
 
 
+def _get_landsea_fraction_fx_dict(variable, config_user):
+    """Get dict of available ``sftlf`` and ``sftof`` variables."""
+    fx_dict = {}
+    fx_vars = ['sftlf']
+    if variable['project'] != 'obs4mips':
+        fx_vars.append('sftof')
+    for fx_var in fx_vars:
+        fx_dict[fx_var] = _get_correct_fx_file(variable, fx_var, config_user)
+    return fx_dict
+
+
 def _update_fx_settings(settings, variable, config_user):
     """Find and set the FX mask settings."""
-    # update for landsea
+    msg = f"Using fx files for %s of dataset {variable['dataset']}:\n%s"
     if 'mask_landsea' in settings:
-        fx_files_dict = {}
-        # Configure ingestion of land/sea masks
         logger.debug('Getting fx mask settings now...')
-        settings['mask_landsea']['fx_files'] = []
-        fx_vars = ['sftlf']
-        if variable['project'] != 'obs4mips':
-            fx_vars.append('sftof')
-        for fx_var in fx_vars:
-            fx_files = _get_correct_fx_file(variable, fx_var, config_user)
-            if fx_files:
-                settings['mask_landsea']['fx_files'].append(fx_files)
+        fx_dict = _get_landsea_fraction_fx_dict(variable, config_user)
+        fx_list = [fx_file for fx_file in fx_dict.values() if fx_file]
+        settings['mask_landsea']['fx_files'] = fx_list
+        logger.info(msg, 'land/sea masking', pformat(fx_dict))
 
     if 'mask_landseaice' in settings:
         logger.debug('Getting fx mask settings now...')
@@ -439,6 +440,13 @@ def _update_fx_settings(settings, variable, config_user):
         if fx_files_dict['sftgif']:
             settings['mask_landseaice']['fx_files'].append(
                 fx_files_dict['sftgif'])
+        logger.info(msg, 'land/sea ice masking', pformat(fx_files_dict))
+
+    if 'weighting_landsea_fraction' in settings:
+        logger.debug("Getting fx files for landsea fraction weighting now...")
+        fx_dict = _get_landsea_fraction_fx_dict(variable, config_user)
+        settings['weighting_landsea_fraction']['fx_files'] = fx_dict
+        logger.info(msg, 'land/sea fraction weighting', pformat(fx_dict))
 
     for step in ('area_statistics', 'volume_statistics'):
         if settings.get(step, {}).get('fx_files'):
@@ -448,6 +456,7 @@ def _update_fx_settings(settings, variable, config_user):
                 fxvar: _get_correct_fx_file(variable, fxvar, config_user)
                 for fxvar in var['fx_files']}
             settings[step]['fx_files'] = fx_files_dict
+            logger.info(msg, step, pformat(fx_files_dict))
 
 
 def _read_attributes(filename):

esmvalcore/preprocessor/__init__.py

@@ -29,6 +29,7 @@
 from ._units import convert_units
 from ._volume import (volume_statistics, depth_integration, extract_trajectory,
                       extract_transect, extract_volume)
+from ._weighting import weighting_landsea_fraction
 
 logger = logging.getLogger(__name__)
 
@@ -53,6 +54,8 @@
     'fix_data',
     # Level extraction
     'extract_levels',
+    # Weighting
+    'weighting_landsea_fraction',
     # Mask landsea (fx or Natural Earth)
     'mask_landsea',
     # Mask landseaice, sftgif only

esmvalcore/preprocessor/_derive/_shared.py

@@ -8,43 +8,6 @@
 logger = logging.getLogger(__name__)
 
 
-def _get_land_fraction(cubes, var_name, derive_from_ocean_fraction=False):
-    """Extract land fraction as :mod:`dask.array`."""
-    cube = cubes.extract_strict(_var_name_constraint(var_name))
-    if derive_from_ocean_fraction:
-        fx_vars = ['sftof', 'sftlf']
-    else:
-        fx_vars = ['sftlf']
-    land_fraction = None
-    for fx_var in fx_vars:
-        if land_fraction is not None:
-            break
-        try:
-            fx_cube = cubes.extract_strict(_var_name_constraint(fx_var))
-        except iris.exceptions.ConstraintMismatchError:
-            logger.debug(
-                "Cannot correct cube '%s' with '%s', fx file not found",
-                var_name, fx_var)
-        else:
-            if not _shape_is_broadcastable(fx_cube.shape, cube.shape):
-                logger.debug("Cannot broadcast fx cube '%s' to cube '%s'",
-                             fx_var, var_name)
-            else:
-                if fx_var == 'sftof':
-                    land_fraction = 1.0 - fx_cube.core_data() / 100.0
-                else:
-                    land_fraction = fx_cube.core_data() / 100.0
-                logger.debug("Using fx cube '%s' to fix '%s'", fx_var,
-                             var_name)
-    return land_fraction
-
-
-def _shape_is_broadcastable(shape_1, shape_2):
-    """Check if two :mod:`numpy.array' shapes are broadcastable."""
-    return all((m == n) or (m == 1) or (n == 1)
-               for (m, n) in zip(shape_1[::-1], shape_2[::-1]))
-
-
 def _var_name_constraint(var_name):
     """:mod:`iris.Constraint` using `var_name` of a :mod:`iris.cube.Cube`."""
     return Constraint(cube_func=lambda c: c.var_name == var_name)
@@ -67,17 +30,3 @@ def cloud_area_fraction(cubes, tau_constraint, plev_constraint):
         new_cube = new_cube.collapsed('air_pressure', iris.analysis.SUM)
 
     return new_cube
-
-
-def grid_area_correction(cubes, var_name, ocean_var=False):
-    """Correct (flux) variable defined relative to land/sea area."""
-    cube = cubes.extract_strict(_var_name_constraint(var_name))
-    core_data = cube.core_data()
-    land_fraction = _get_land_fraction(cubes,
-                                       var_name,
-                                       derive_from_ocean_fraction=ocean_var)
-    if land_fraction is not None:
-        if ocean_var:
-            land_fraction = 1.0 - land_fraction
-        cube.data = core_data * land_fraction
-    return cube