docs/notebooks/units.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unit handling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "\n",
    "import xclim as xc\n",
    "\n",
    "# Set display to HTML style (optional)\n",
    "xr.set_options(display_style=\"html\", display_width=50)\n",
    "\n",
    "# import plotting stuff\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "plt.style.use(\"seaborn\")\n",
    "plt.rcParams[\"figure.figsize\"] = (11, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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 metadata attribute for units\n",
    "\n",
    "For precipitation data, `xclim` expects precipitation fluxes. This could be units of `length/time`, such as mm/d, or units of `mass / area / time`, such as kg/m²/s. Units of `length` only, such as mm, are not supported, because the interpretation depends on the frequency of the data, which cannot always be inferred explicitly from the data. For example, if a daily precipitation series records total daily precipitation and has units of `mm`, change the units attribute to `mm/d` before computing indicators. Note that `xclim` will automatically convert between `mass / area / time` and `length/time` using a water density of 1000 kg/m³ when the context is hydrology.\n",
    "\n",
    "In the following examples, our toy temperature dataset comes in units of Kelvins (`\"degK\"`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# See the Usage page for details on opening datasets, subsetting and resampling.\n",
    "ds = xr.tutorial.open_dataset(\"air_temperature\")\n",
    "tas = (\n",
    "    ds.air.sel(lat=40, lon=270, method=\"nearest\")\n",
    "    .resample(time=\"D\")\n",
    "    .mean(keep_attrs=True)\n",
    ")\n",
    "print(tas.attrs[\"units\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using [pint](https://pint.readthedocs.io/), `xclim` provides useful functions to convert the units of datasets and `DataArray`s. Here, we convert our kelvin data to the very useful Fahrenheits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tas_F = xc.units.convert_units_to(tas, \"degF\")\n",
    "print(tas_F.attrs[\"units\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Threshold indices\n",
    "\n",
    "`xclim` unit handling also applies to threshold indicators. Users can provide threshold in units of choice and `xclim` will adjust automatically. For example determining 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. Alernatively users can even provide a threshold in Kelvin '293.15 K' (if they really wanted to)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with xc.set_options(cf_compliance=\"log\"):\n",
    "    # Using Kelvin data, threshold in Celsius\n",
    "    out1 = xc.atmos.tx_days_above(tasmax=tas, thresh=\"20 C\", freq=\"MS\")\n",
    "\n",
    "    # Using Fahrenheit data, threshold in Celsius\n",
    "    out2 = xc.atmos.tx_days_above(tasmax=tas_F, thresh=\"20 C\", freq=\"MS\")\n",
    "\n",
    "    # Using Fahrenheit data, with threshold in Kelvin\n",
    "    out3 = xc.atmos.tx_days_above(tasmax=tas_F, thresh=\"293.15 K\", freq=\"MS\")\n",
    "\n",
    "# Plot and see that it's all identical:\n",
    "plt.figure()\n",
    "out1.plot(label=\"K and degC\", linestyle=\"-\")\n",
    "out2.plot(label=\"degF and degC\", marker=\"s\", markersize=10, linestyle=\"none\")\n",
    "out3.plot(label=\"degF and K\", marker=\"o\", linestyle=\"none\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sum and count indices\n",
    "\n",
    "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.\n",
    "\n",
    "Indicators, on the other hand, do not have this flexibility and often **expect** input at a given frequency, more often daily then otherwise.\n",
    "\n",
    "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)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tas_6h = ds.air.sel(\n",
    "    lat=40, lon=270, method=\"nearest\"\n",
    ")  # no resampling, original data is 6-hourly\n",
    "out4_h = xc.indices.tx_days_above(tasmax=tas_6h, thresh=\"20 C\", freq=\"MS\")\n",
    "out4_h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out4_d = xc.units.convert_units_to(out4_h, \"d\")\n",
    "plt.figure()\n",
    "out1.plot(label=\"From daily input\", linestyle=\"-\")\n",
    "out4_d.plot(label=\"From 6-hourly input\", linestyle=\"-\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other utilites\n",
    "\n",
    "Many helper functions are defined in `xclim.core.units`, see [Unit handling module](../api.rst#module-xclim.core.units)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}