Unit Handling

[ ]:

from __future__ import annotations

import xarray as xr

import xclim
from xclim import indices
from xclim.core import units
from xclim.testing import open_dataset

# Set display to HTML style (optional)
xr.set_options(display_style="html", display_width=50)

# import plotting stuff
import matplotlib.pyplot as plt
import nc_time_axis  # noqa

%matplotlib inline
plt.rcParams["figure.figsize"] = (11, 5)

A lot of effort has been placed into automatic handling of input data units. xclim will automatically detect the input variable(s) units (e.g. °C versus K or mm/s versus mm/day etc.) and adjust on-the-fly in order to calculate indices in the consistent manner. This comes with the obvious caveat that input data requires a metadata attribute for units : the units attribute is required, and the standard_name can be useful for automatic conversions.

The main unit handling method is `xclim.core.units.convert_units_to <../xclim.core.rst#xclim.core.units.convert_units_to>`__ which can also be useful on its own. xclim relies on pint for unit handling and extends the units registry and formatting functions of cf-xarray.

Simple example: Temperature

[ ]:

# See the Usage page for details on opening datasets, subsetting and resampling.
ds = xr.tutorial.open_dataset("air_temperature")
tas = (
    ds.air.sel(lat=40, lon=270, method="nearest")
    .resample(time="D")
    .mean(keep_attrs=True)
)
print(tas.attrs["units"])

Here, we convert our kelvin data to the very useful Fahrenheits:

[ ]:

tas_F = units.convert_units_to(tas, "degF")
print(tas_F.attrs["units"])

Smart conversions: Precipitation

For precipitation data, xclim usually expects precipitation fluxes, so units of mass / (area * time). However, many indicators will also accept rates (length / time, for example mm/d) by implicitly assuming the data refers to liquid water, and thus that we can simply multiply by 1000 kg/m³ to convert from the latter to the former. This transformation is enabled on indicators that have Indicator.context == 'hydro'.

We can also leverage the CF-conventions to perform some other “smart” conversions. For example, if the CF standard name of an input refers to liquid water, the flux ⇋ rate and amount ⇋ thickness conversions explained above will be automatic in xc.core.units.convert_units_to, whether the “hydro” context is activated or not. Another CF-driven conversion is between amount and flux or thickness and rate. Here again, convert_units_to will see if the standard_name attribute, but it will also need to infer the frequency of the data. For example, if a daily precipitation series records total daily precipitation and has units of mm (a “thickness”), it should use the lwe_thickness_of_precipitation_amount standard name and have a regular time coordinate, With these two, xclim will understand it and be able to convert it to a precipitation flux (by dividing by 1 day and multiplying by 1000 kg/m³).

These CF conversions are not automatically done in the indicator (in opposition to the “hydro” context ones). convert_units_to should be called beforehand.

Here are some examples:

Going from a precipitation flux to a daily thickness

[ ]:

ds = open_dataset("ERA5/daily_surface_cancities_1990-1993.nc")
ds.pr.attrs

[ ]:

pr_as_daily_total = units.convert_units_to(ds.pr, "mm")
pr_as_daily_total

Going from the liquid water equivalent thickness of snow (swe) to the snow amount (snw).

The former being common in observations and the latter being the CMIP6 variable. Notice that convert_units_to cannot handle the variable name itself, that change has to be done by hand.

[ ]:

ds.swe.attrs

[ ]:

units.convert_units_to(ds.swe, "kg m-2").rename("snw")

Threshold indices

xclim unit handling also applies to threshold indicators. Users can provide threshold in units of choice and xclim will adjust automatically. For example, in order to determine the number of days with tasmax > 20 °C, users can define a threshold input of "20 C" or "20 degC" even if input data is in Kelvin. Alternatively, users can even provide a threshold in Kelvin ("293.15 K"), if they really wanted to.

[ ]:

tasmaxK = ds.tasmax.sel(location="Halifax")
tasmaxF = units.convert_units_to(ds.tasmax.sel(location="Halifax"), "degF")

with xclim.set_options(cf_compliance="log"):
    # Using Kelvin data, threshold in Celsius
    out1 = xclim.atmos.tx_days_above(tasmax=tasmaxK, thresh="20 C", freq="YS")

    # Using Fahrenheit data, threshold in Celsius
    out2 = xclim.atmos.tx_days_above(tasmax=tasmaxF, thresh="20 C", freq="YS")

    # Using Fahrenheit data, with threshold in Kelvin
    out3 = xclim.atmos.tx_days_above(tasmax=tasmaxF, thresh="293.15 K", freq="YS")

# Plot and see that it's all identical:
plt.figure()
out1.plot(label="K and degC", linestyle="-")
out2.plot(label="degF and degC", marker="s", markersize=10, linestyle="none")
out3.plot(label="degF and K", marker="o", linestyle="none")
plt.legend()

Sum and count indices

Many indices in xclim will either sum values or count events along the time dimension and over a period. As of version 0.24, unit handling dynamically infers what the sampling frequency and its corresponding unit is.

Indicators, on the other hand, do not have this flexibility and often expect input at a given frequency, more often daily than otherwise.

For example, we can run the tx_days_above on the 6-hourly test data, and it should return similar results as on the daily data, but in units of h (the base unit of the sampling frequency).

[ ]:

ds = xr.tutorial.open_dataset("air_temperature")
tas_6h = ds.air.sel(
    lat=40, lon=270, method="nearest"
)  # no resampling, original data is 6-hourly
tas_D = tas_6h.resample(time="D").mean()
out1_h = xclim.indices.tx_days_above(tasmax=tas_6h, thresh="20 C", freq="MS")
out2_D = xclim.indices.tx_days_above(tasmax=tas_D, thresh="20 C", freq="MS")
out1_h

[ ]:

out1_D = units.convert_units_to(out1_h, "d")
plt.figure()
out2_D.plot(label="From daily input", linestyle="-")
out1_D.plot(label="From 6-hourly input", linestyle="-")
plt.legend()

Other utilities

Many helper functions are defined in xclim.core.units, see Unit handling module.