""" Measures Submodule ================== Measures compare adjusted simulations to a reference SDBA diagnostic tests are made up of properties and measures. Measures compare adjusted simulations to a reference, through statistical properties or directly. This framework for the diagnostic tests was inspired by the `VALUE `_ project. This module depends on `xclim`. Run `pip install xsdba['extras']` to install it. """ from __future__ import annotations from collections.abc import Sequence import numpy as np import xarray as xr from xclim.core.indicator import Indicator, base_registry from xsdba.base import Grouper from xsdba.typing import InputKind from xsdba.units import convert_units_to, pint2cfattrs, units2pint from xsdba.utils import _pairwise_spearman class StatisticalMeasure(Indicator): """ Base indicator class for statistical measures used when validating bias-adjusted outputs. Statistical measures use input data where the time dimension was reduced, usually by the computation of a :py:class:`xsdba.properties.StatisticalProperty` instance. They usually take two arrays as input: "sim" and "ref", "sim" being measured against "ref". The two arrays must have identical coordinates on their common dimensions. Statistical measures are generally unit-generic. If the inputs have different units, "sim" is converted to match "ref". """ realm = "generic" @classmethod def _ensure_correct_parameters(cls, parameters): inputs = {k for k, p in parameters.items() if p.kind == InputKind.VARIABLE} if not inputs.issuperset({"sim", "ref"}): raise ValueError(f"{cls.__name__} requires 'sim' and 'ref' as inputs. Got {inputs}.") return super()._ensure_correct_parameters(parameters) def _preprocess_and_checks(self, das, params): """Perform parent's checks and also check convert units so that sim matches ref.""" das, params = super()._preprocess_and_checks(das, params) # Convert grouping and check if allowed: das["sim"] = convert_units_to(das["sim"], das["ref"]) ref = das["ref"] sim = das["sim"] # Check if common coordinates are identical. newsim, newref = xr.broadcast(sim, ref) for dim in set(sim.dims).union(ref.dims): if [sim[dim].size, ref[dim].size] != [newsim[dim].size, newref[dim].size]: raise ValueError(f"Common dimension {dim} has different coordinates between ref and sim.") return das, params class StatisticalPropertyMeasure(Indicator): """ Base indicator class for statistical properties that include the comparison measure, used when validating bias-adjusted outputs. StatisticalPropertyMeasure objects combine the functionalities of :py:class:`xsdba.properties.StatisticalProperty` and :py:class:`xsdba.properties.StatisticalMeasure`. Statistical properties usually reduce the time dimension and sometimes more dimensions (for example in spatial properties), sometimes adding a grouping dimension according to the passed value of `group` (e.g.: group='time.month' means the loss of the time dimension and the addition of a month one). Statistical measures usually take two arrays as input: "sim" and "ref", "sim" being measured against "ref". Statistical property-measures are generally unit-generic. If the inputs have different units, "sim" is converted to match "ref". """ aspect = None """The aspect the statistical property studies: marginal, temporal, multivariate or spatial.""" allowed_groups = None """A list of allowed groupings. A subset of dayofyear, week, month, season or group. The latter stands for no temporal grouping.""" realm = "generic" @classmethod def _ensure_correct_parameters(cls, parameters): inputs = {k for k, p in parameters.items() if p.kind == InputKind.VARIABLE} if not inputs.issuperset({"sim", "ref"}): raise ValueError(f"{cls.__name__} requires 'sim' and 'ref' as inputs. Got {inputs}.") if "group" not in parameters: raise ValueError( f"{cls.__name__} require a 'group' argument, use the base Indicator class if your computation doesn't perform any regrouping." ) return super()._ensure_correct_parameters(parameters) def _preprocess_and_checks(self, das, params): """Perform parent's checks and also check convert units so that sim matches ref.""" das, params = super()._preprocess_and_checks(das, params) das["sim"] = convert_units_to(das["sim"], das["ref"]) # Convert grouping and check if allowed: if isinstance(params["group"], str): params["group"] = Grouper(params["group"]) if self.allowed_groups is not None: if params["group"].prop not in self.allowed_groups: raise ValueError( f"Grouping period {params['group'].prop_name} is not allowed for property " f"{self.identifier} (needs something in " f"{list(map(lambda g: '.' + g.replace('group', ''), self.allowed_groups))})." ) return das, params def _postprocess(self, outs, das, params): """Squeeze `group` dim if needed.""" outs = super()._postprocess(outs, das, params) for ii, out in enumerate(outs): if "group" in out.dims: outs[ii] = out.squeeze("group", drop=True) return outs base_registry["StatisticalMeasure"] = StatisticalMeasure base_registry["StatisticalPropertyMeasure"] = StatisticalPropertyMeasure def _bias(sim: xr.DataArray, ref: xr.DataArray) -> xr.DataArray: """ Bias. The bias is the simulation minus the reference. Parameters ---------- sim : xr.DataArray data from the simulation (one value for each grid-point). ref : xr.DataArray data from the reference (observations) (one value for each grid-point). Returns ------- xr.DataArray, [same as ref] Absolute bias. """ out = sim - ref out.attrs.update(pint2cfattrs(units2pint(ref.attrs["units"]), is_difference=True)) return out bias = StatisticalMeasure(identifier="bias", compute=_bias) def _relative_bias(sim: xr.DataArray, ref: xr.DataArray) -> xr.DataArray: """ Relative Bias. The relative bias is the simulation minus reference, divided by the reference. Parameters ---------- sim : xr.DataArray data from the simulation (one value for each grid-point). ref : xr.DataArray data from the reference (observations) (one value for each grid-point). Returns ------- xr.DataArray, [dimensionless] Relative bias. """ out = (sim - ref) / ref return out.assign_attrs(units="") relative_bias = StatisticalMeasure(identifier="relative_bias", compute=_relative_bias, units="") def _circular_bias(sim: xr.DataArray, ref: xr.DataArray) -> xr.DataArray: """ Circular bias. Bias considering circular time series. E.g. The bias between doy 365 and doy 1 is 364, but the circular bias is -1. Parameters ---------- sim : xr.DataArray data from the simulation (one value for each grid-point). ref : xr.DataArray data from the reference (observations) (one value for each grid-point). Returns ------- xr.DataArray, [days] Circular bias. """ out = (sim - ref) % 365 out = out.where(out <= 365 / 2, 365 - out) # when condition false, replace by 2nd arg out = out.where(ref >= sim, out * -1) # when condition false, replace by 2nd arg return out.assign_attrs(units="days") circular_bias = StatisticalMeasure(identifier="circular_bias", compute=_circular_bias, units="days") def _ratio(sim: xr.DataArray, ref: xr.DataArray) -> xr.DataArray: """ Ratio. The ratio is the quotient of the simulation over the reference. Parameters ---------- sim : xr.DataArray data from the simulation (one value for each grid-point). ref : xr.DataArray data from the reference (observations) (one value for each grid-point). Returns ------- xr.DataArray, [dimensionless] Ratio. """ out = sim / ref out.attrs["units"] = "" return out ratio = StatisticalMeasure(identifier="ratio", compute=_ratio, units="") def _rmse(sim: xr.DataArray, ref: xr.DataArray, group: str | Grouper = "time") -> xr.DataArray: """ Root mean square error. The root mean square error on the time dimension between the simulation and the reference. Parameters ---------- sim : xr.DataArray Data from the simulation (a time-series for each grid-point). ref : xr.DataArray Data from the reference (observations) (a time-series for each grid-point). group: str Compute the property and measure for each temporal groups individually. Currently not implemented. Returns ------- xr.DataArray, [same as ref] Root mean square error. """ def _rmse_internal(_sim: xr.DataArray, _ref: xr.DataArray) -> xr.DataArray: _f: xr.DataArray = np.sqrt(np.nanmean((_sim - _ref) ** 2, axis=-1)) return _f out = xr.apply_ufunc( _rmse_internal, sim, ref, input_core_dims=[["time"], ["time"]], dask="parallelized", ) out = out.assign_attrs(pint2cfattrs(units2pint(ref.units), is_difference=True)) return out rmse = StatisticalPropertyMeasure( identifier="rmse", aspect="temporal", compute=_rmse, allowed_groups=["group"], cell_methods="time: mean", ) def _mae(sim: xr.DataArray, ref: xr.DataArray, group: str | Grouper = "time") -> xr.DataArray: """ Mean absolute error. The mean absolute error on the time dimension between the simulation and the reference. Parameters ---------- sim : xr.DataArray data from the simulation (a time-series for each grid-point). ref : xr.DataArray data from the reference (observations) (a time-series for each grid-point). group : str Compute the property and measure for each temporal groups individually. Currently not implemented. Returns ------- xr.DataArray, [same as ref] Mean absolute error. """ def _mae_internal(_sim: xr.DataArray, _ref: xr.DataArray) -> xr.DataArray: _f: xr.DataArray = np.nanmean(np.abs(_sim - _ref), axis=-1) return _f out = xr.apply_ufunc( _mae_internal, sim, ref, input_core_dims=[["time"], ["time"]], dask="parallelized", ) out = out.assign_attrs(pint2cfattrs(units2pint(ref.units), is_difference=True)) return out mae = StatisticalPropertyMeasure( identifier="mae", aspect="temporal", compute=_mae, allowed_groups=["group"], cell_methods="time: mean", ) def _annual_cycle_correlation( sim: xr.DataArray, ref: xr.DataArray, window: int = 15, group: str | Grouper = "time", ) -> xr.DataArray: """ Annual cycle correlation. Pearson correlation coefficient between the smooth day-of-year averaged annual cycles of the simulation and the reference. In the smooth day-of-year averaged annual cycles, each day-of-year is averaged over all years and over a window of days around that day. Parameters ---------- sim : xr.DataArray data from the simulation (a time-series for each grid-point). ref : xr.DataArray data from the reference (observations) (a time-series for each grid-point). window : int Size of window around each day of year around which to take the mean. E.g. If window=31, Jan 1st is averaged over from December 17th to January 16th. group : str Compute the property and measure for each temporal groups individually. Currently not implemented. Returns ------- xr.DataArray, [dimensionless] Annual cycle correlation. """ # group by day-of-year and window around each doy grouper_test = Grouper("time.dayofyear", window=window) # for each day, mean over X day window and over all years to create a smooth avg annual cycle sim_annual_cycle = grouper_test.apply("mean", sim) ref_annual_cycle = grouper_test.apply("mean", ref) out = xr.corr(ref_annual_cycle, sim_annual_cycle, dim="dayofyear") return out.assign_attrs(units="") annual_cycle_correlation = StatisticalPropertyMeasure( identifier="annual_cycle_correlation", aspect="temporal", compute=_annual_cycle_correlation, allowed_groups=["group"], ) def _scorr( sim: xr.DataArray, ref: xr.DataArray, *, dims: Sequence | None = None, group: str | Grouper = "time", ): """ Spatial correllogram. Compute the inter-site correlations of each array, compute the difference in correlations and sum. Taken from Vrac (2018). The spatial and temporal dimensions are reduced. Parameters ---------- sim : xr.DataArray data from the simulation (a time-series for each grid-point). ref : xr.DataArray data from the reference (observations) (a time-series for each grid-point). dims : sequence of strings, optional Name of the spatial dimensions. If None (default), all dimensions except 'time' are used. group : str Compute the property and measure for each temporal groups individually. Currently not implemented. Returns ------- xr.DataArray, [dimensionless] Sum of the inter-site correlation differences. """ if dims is None: dims = [d for d in sim.dims if d != "time"] refcorr = _pairwise_spearman(ref, dims) simcorr = _pairwise_spearman(sim, dims) S_corr = (simcorr - refcorr).sum(["_spatial", "_spatial2"]) return S_corr.assign_attrs(units="") scorr = StatisticalPropertyMeasure(identifier="Scorr", aspect="spatial", compute=_scorr, allowed_groups=["group"]) def _taylordiagram( sim: xr.DataArray, ref: xr.DataArray, dim: str = "time", group: str | Grouper = "time", normalize: bool = False, ) -> xr.DataArray: """ Taylor diagram. Compute the respective standard deviations of a simulation and a reference array, as well as the Pearson correlation coefficient between both, all necessary parameters to plot points on a Taylor diagram. Parameters ---------- sim : xr.DataArray data from the simulation (a time-series for each grid-point). ref : xr.DataArray data from the reference (observations) (a time-series for each grid-point). dim : str Dimension across which the correlation and standard deviation should be computed. group : str Compute the property and measure for each temporal groups individually. Currently not implemented. normalize : bool If `True`, divide the standard deviations by the standard deviation of the reference. Default is `False`. Returns ------- xr.DataArray, [same as ref] Standard deviations of sim, ref and correlation coefficient between both. """ corr = xr.corr(sim, ref, dim=dim) ref_std = ref.std(dim=dim, skipna=True, keep_attrs=True) sim_std = sim.std(dim=dim, skipna=True, keep_attrs=True) new_dim = xr.DataArray(["ref_std", "sim_std", "corr"], dims=("taylor_param",), name="taylor_param") out = xr.concat( [ref_std, sim_std, corr], new_dim, coords="minimal", compat="override", # Take common coords from `ref_std`. ).assign_attrs( { "correlation_type": "Pearson correlation coefficient", "units": ref.units, } ) # Normalize the standard deviations byt the standard deviation of the reference. if normalize: if (out[{"taylor_param": 0}] == 0).any(): raise ValueError("`ref_std =0` (homogeneous field) obtained, normalization is not possible.") with xr.set_options(keep_attrs=True): out[{"taylor_param": [0, 1]}] = out[{"taylor_param": [0, 1]}] / out[{"taylor_param": 0}] out.attrs["normalized"] = True out.attrs["units"] = "" return out taylordiagram = StatisticalPropertyMeasure( identifier="taylordiagram", aspect="temporal", compute=_taylordiagram, allowed_groups=["group"], )