docs/notebooks/example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "keep_output": false
   },
   "source": [
    "# Workflow Examples\n",
    "\n",
    "\n",
    "`xclim` is built on very powerful multiprocessing and distributed computation libraries, notably `xarray` and `dask`.\n",
    "\n",
    "`xarray` is a python package making it easy to work with n-dimensional arrays. It labels axes with their names `[time, lat, lon, level]` instead of indices `[0,1,2,3]`, reducing the likelihood of bugs and making the code easier to understand. One of the key strengths of `xarray` is that it knows how to deal with non-standard calendars (we're looking at you, \"`360_days`\") and can easily resample daily time series to weekly, monthly, seasonal or annual periods.  Finally, `xarray` is tightly integrated with `dask`, a package that can automatically parallelize operations.\n",
    "\n",
    "The following are a few examples to consult when using `xclim` to subset NetCDF arrays and compute climate indicators, taking advantage of the parallel processing capabilities offered by `xarray` and `dask`. For more information about these projects, please see their documentation pages:\n",
    "\n",
    "* [xarray documentation](https://xarray.pydata.org/en/stable/)\n",
    "* [dask documentation](https://docs.dask.org/en/stable/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "keep_output": false
   },
   "source": [
    "## Environment configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "from __future__ import annotations\n",
    "\n",
    "# File handling libraries\n",
    "import tempfile\n",
    "from pathlib import Path\n",
    "\n",
    "import xarray as xr\n",
    "\n",
    "# Set xarray to use HTML for displaying outputs\n",
    "xr.set_options(display_style=\"html\")\n",
    "\n",
    "import xclim.indices\n",
    "\n",
    "# Output folder\n",
    "output_folder = Path(tempfile.mkdtemp())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "keep_output": false
   },
   "source": [
    "## Setting up the Dask client: parallel processing\n",
    "\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "In this example, we are using the ``**dask.distributed**`` submodule. This is not installed by default in a basic `xclim` installation. Be sure to add `distributed` to your Python installation before setting up parallel processing operations!\n",
    "\n",
    "</div>\n",
    "\n",
    "First, we create a pool of workers that will wait for jobs. The `xarray` library will automatically connect to these workers and dispatch them jobs that can be run in parallel.\n",
    "\n",
    "The dashboard link lets you see in real time how busy those workers are.\n",
    "\n",
    "* [dask distributed documentation](https://distributed.dask.org/en/latest/)\n",
    "\n",
    "This step is not mandatory, as `dask` will fall back to its \"single machine scheduler\" if a `Client` is not created. However, this default scheduler doesn't allow you to set the number of threads or a memory limit and doesn't start the dashboard, which can be quite useful to understand your task's progress.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "from distributed import Client\n",
    "\n",
    "# Depending on your workstation specifications, you may need to adjust these values.\n",
    "# On a single machine, n_workers=1 is usually better.\n",
    "client = Client(n_workers=1, threads_per_worker=4, memory_limit=\"4GB\")\n",
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "keep_output": false
   },
   "source": [
    "## Creating xarray datasets\n",
    "\n",
    "To open a NetCDF file with `xarray`, we use `xr.open_dataset(<path to file>)`. By default, the entire file is stored in one chunk, so there is no parallelism. To trigger parallel computations, we need to explicitly specify the **chunk size**.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "In this example, instead of opening a local file, we pass an *OPeNDAP* URL to xarray. It retrieves the data automatically. Notice also that opening the dataset is quite fast. In fact, the data itself has not been downloaded yet, only the coordinates and the metadata. The downloads will be triggered only when the values need to be accessed directly.\n",
    "\n",
    "</div>\n",
    "\n",
    "`dask`'s parallelism is based on memory chunks; We need to tell `xarray` to split our NetCDF array into chunks of a given size, and operations on each chunk of the array will automatically be dispatched to the workers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "data_url = \"https://pavics.ouranos.ca/twitcher/ows/proxy/thredds/dodsC/datasets/simulations/bias_adjusted/cmip5/ouranos/cb-oura-1.0/day_ACCESS1-3_historical+rcp85_r1i1p1_na10kgrid_qm-moving-50bins-detrend_1950-2100.ncml\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "# Chunking in memory along the time dimension.\n",
    "# Note that the data type is a 'dask.array'. xarray will automatically use client workers.\n",
    "ds = xr.open_dataset(\n",
    "    data_url,\n",
    "    chunks={\"time\": 365, \"lat\": 168, \"lon\": 150},\n",
    "    drop_variables=[\"ts\", \"time_vectors\"],\n",
    ")\n",
    "print(ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "print(ds.tasmin.chunks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "keep_output": false
   },
   "source": [
    "## Multi-file datasets\n",
    "\n",
    "NetCDF files are often split into periods to keep file size manageable. A single dataset can be split in dozens of individual files. `xarray` has a function `open_mfdataset` that can open and aggregate a list of files and construct a unique *logical* dataset. `open_mfdataset` can aggregate files over coordinates (time, lat, lon) and variables.\n",
    "\n",
    "* Note that opening a multi-file dataset automatically chunks the array (one chunk per file).\n",
    "* Note also that because `xarray` reads every file metadata to place it in a logical order, it can take a while to load."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "## Create multi-file data & chunks\n",
    "# ds = xr.open_mfdataset('/path/to/files*.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Subsetting and selecting data with xarray\n",
    "Here, we will reduce the size of our data using the methods implemented in xarray ([docs here](https://docs.xarray.dev/en/stable/user-guide/indexing.html))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "ds2 = ds.sel(lat=slice(50, 45), lon=slice(-70, -65), time=slice(\"2090\", \"2100\"))\n",
    "print(ds2.tasmin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds3 = ds.sel(lat=46.8, lon=-71.22, method=\"nearest\").sel(time=\"1993\")\n",
    "print(ds3.tasmin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more powerful subsetting tools with features such as coordinate reference system (CRS) aware subsetting and vector shape masking, the `xclim` developers strongly encourage users to consider the subsetting utilities of the [clisops](https://github.com/roocs/clisops) package.\n",
    "\n",
    "Their documentation showcases several examples of how to perform more complex subsetting: [clisops.core.subset](https://clisops.readthedocs.io/en/latest/notebooks/core_subset.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "keep_output": false
   },
   "source": [
    "## Climate index calculation & resampling frequencies\n",
    "\n",
    "`xclim` has two layers for the calculation of indicators. The bottom layer is composed of a list of functions that take one or more `xarray.DataArray`'s as input and return an `xarray.DataArray` as output. You'll find these functions in `xclim.indices`. The indicator's logic is contained in this function, as well as some unit handling, but it doesn't perform any data consistency checks (like if the time frequency is daily), and doesn't adjust the metadata of the output array.\n",
    "\n",
    "The second layer are class instances that you'll find organized by *realm*. So far, there are three realms available in `xclim.atmos`, `xclim.seaIce` and `xclim.land`, the first one being the most exhaustive. Before running computations, these classes check if the input data is a daily average of the expected variable:\n",
    "\n",
    "1. If an indicator expects a daily mean, and you pass it a daily max, a `warning` will be raised.\n",
    "2. After the computation, it also checks the number of values per period to make sure there are not missing values or `NaN` in the input data. If there are, the output is going to be set to `NaN`.  Ex. : If the indicator performs a yearly resampling, but there are only 350 non-`NaN` values in one given year in the input data, that year's output will be `NaN`.\n",
    "3. The output units are set correctly as well as other properties of the output array, complying as much as possible with CF conventions.\n",
    "\n",
    "For new users, we suggest you use the classes found in `xclim.atmos` and others. If you know what you're doing, and you want to circumvent the built-in checks, then you can use the `xclim.indices` directly.\n",
    "\n",
    "Almost all `xclim` indicators convert daily data to lower time frequencies, such as seasonal or annual values. This is done using `xarray.DataArray.resample` method. Resampling creates a grouped object over which you apply a reduction operation (e.g. mean, min, max). The list of available frequency is given in the link below, but the most often used are:\n",
    "\n",
    "- `YS`: annual starting in January\n",
    "- `YS-JUL`: annual starting in July\n",
    "- `MS`: monthly\n",
    "- `QS-DEC`: seasonal starting in December\n",
    "\n",
    "More info about this specification can be found in [pandas' documentation](http://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html#timeseries-offset-aliases)\n",
    "\n",
    "Note - not all offsets in the link are supported by ``cftime`` objects in `xarray`.\n",
    "\n",
    "\n",
    "In the example below, we're computing the **annual maximum temperature of the daily maximum temperature (tx_max)**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "out = xclim.atmos.tx_max(ds2.tasmax, freq=\"YS\")\n",
    "print(out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "keep_output": false
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "If you execute the cell above, you'll see that this operation is quite fast. This a feature coming from dask. Read *Lazy computation*  further down.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "keep_output": false
   },
   "source": [
    "### Comparison of **atmos** vs **indices** modules\n",
    "Using the `xclim.indices` module performs not checks and only fills the `units` attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "out = xclim.indices.tx_days_above(ds2.tasmax, thresh=\"30 degC\", freq=\"YS\")\n",
    "print(out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "keep_output": false
   },
   "source": [
    "With `xclim.atmos`, checks are performed and many CF-compliant attributes are added:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "out = xclim.atmos.tx_days_above(ds2.tasmax, thresh=\"30 degC\", freq=\"YS\")\n",
    "print(out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "# We have created an xarray data-array - We can insert this into an output xr.Dataset object with a copy of the original dataset global attrs\n",
    "ds_out = xr.Dataset(attrs=ds2.attrs)\n",
    "\n",
    "# Add our climate index as a data variable to the dataset\n",
    "ds_out[out.name] = out\n",
    "print(ds_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Different ways of resampling\n",
    "\n",
    "Many indices use algorithms that find the length of given sequences. For instance, `xclim.indices.heat_wave_max_length` finds the longest sequence where `tasmax` and `tasmin` are above given threshold values. Resampling can be used to find the longest sequence in given periods of time, for instance the longest heat wave for each month if the resampling frequency is `freq == \"MS\"`.\n",
    "\n",
    "The order of the two operations just described, i.e. :\n",
    "* Finding the length of sequences respecting a certain criterion (\"run length algorithms\")\n",
    "* Separating the dataset in given time periods (\"resampling\")\n",
    "\n",
    "is important and can lead to differing results.\n",
    "\n",
    "The cell below illustrates this by looking at the maximum lengths of heat waves each month from May 2010 to August 2010 by doing these operations in the two possible orders. The heat wave max lengths for July in a small region of interest $\\text{lat} \\in [43, 44.5],\\, \\text{lon} \\in [-117.5, -116]$ are shown: The maximal lengths are sometimes longer first applying the run length algorithm (`resample_before_rl == False`).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "# import plotting stuff\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (11, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "ds_reduced = (\n",
    "    ds.sel(lat=slice(44.5, 43))\n",
    "    .sel(lon=slice(-117.5, -116))\n",
    "    .sel(time=slice(\"2010-05-01\", \"2011-08-31\"))\n",
    ")\n",
    "tn, tx = ds_reduced.tasmin, ds_reduced.tasmax\n",
    "freq = \"MS\"\n",
    "\n",
    "# Threshold on tasmin: chosen smaller than the default 22.0 degC\n",
    "thresh_tn = \"20.0 degC\"\n",
    "\n",
    "# Computing index by resampling **before** run length algorithm (default value)\n",
    "hw_before = xclim.indices.heat_wave_max_length(\n",
    "    tn, tx, freq=freq, thresh_tasmin=thresh_tn, resample_before_rl=True\n",
    ")\n",
    "# Computing index by resampling **after** run length algorithm\n",
    "hw_after = xclim.indices.heat_wave_max_length(\n",
    "    tn, tx, freq=freq, thresh_tasmin=thresh_tn, resample_before_rl=False\n",
    ")\n",
    "\n",
    "hw_before.sel(time=\"2010-07-01\").plot(vmin=0, vmax=7)\n",
    "plt.title(\"Resample, then run length\")\n",
    "plt.figure()\n",
    "hw_after.sel(time=\"2010-07-01\").plot(vmin=0, vmax=7)\n",
    "plt.title(\"Run length, then resample\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's focus on the point $(-117.21075, 44.29087)$,  which has a maximum wave length of four (4) or seven (7), depending on whether resampling occurs _before_ or _after_ the run length algorithm.\n",
    "\n",
    "Plotting the values of `tasmin` in July and early August, we see a sequence of seven hot minimal temperatures at the end of July that surpass the threshold to qualify for a heat wave.\n",
    "\n",
    "1. If resampling occurs first, and we first separate the periods in months, the run length algorithms will only look for sequences of hot days **within** the month of July and will exclude the last 3 days of this sequence of 7 days.\n",
    "\n",
    "2. Using the run length algorithm before resampling looks for sequences of hot days in all the dataset given (temperatures from May 1, 2010 to Aug. 31, 2010) and then subdivides these sequences in the months where they have started. Since it starts in July, this sequence is registered as counts for a heat wave of seven days happening in July.\n",
    "\n",
    "This also implies that the first 3 days of August which belong in this sequence of seven days will be counted as a sequence in August with the first method, but not with the second."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "from xclim.core.units import convert_units_to\n",
    "\n",
    "# Select a spatial point of interest in July-early August\n",
    "lon_i, lat_i = -117.21075, 44.29087\n",
    "tn_pt = tn.sel(time=slice(\"2010-07-01\", \"2010-08-05\")).sel(lat=lat_i, lon=lon_i)\n",
    "tn_pt = convert_units_to(tn_pt, \"degC\")\n",
    "\n",
    "tn_pt.plot(marker=\"o\", label=\"tasmin\")\n",
    "plt.axhline(y=convert_units_to(thresh_tn, \"degC\"), color=\"orange\", label=thresh_tn)\n",
    "plt.axvline(x=[\"2010-08-01\"], color=\"green\", label=\"Aug. 1st\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "keep_output": false
   },
   "source": [
    "## *Lazy* computation - Nothing has been computed so far !\n",
    "\n",
    "If you look at the output of those operations, they're identified as `dask.array` objects. What happens is that `dask` creates a chain of operations that, when executed, will yield the values we want. We have thus far only created a schedule of tasks with a small preview and not done any actual computations. You can trigger computations by using the `load` or `compute` method, or writing the output to disk via `to_netcdf`. Of course, calling `.plot()` will also trigger the computation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "output_file = output_folder / \"test_tx_max.nc\"\n",
    "dsOut.to_netcdf(output_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "keep_output": false
   },
   "source": [
    "*(Times may of course vary depending on the machine and the Client settings)*\n",
    "\n",
    "### Performance tips\n",
    "#### Optimizing the chunk size\n",
    "\n",
    "You can improve performance by being smart about chunk sizes. If chunks are too small, there is a lot of time lost in overhead. If chunks are too large, you may end up exceeding the individual worker memory limit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "print(ds2.chunks[\"time\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "# rechunk data in memory for the entire grid\n",
    "ds2c = ds2.chunk(chunks={\"time\": 4 * 365})\n",
    "print(ds2c.chunks[\"time\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "out = xclim.atmos.tx_max(ds2c.tasmax, freq=\"YS\")\n",
    "ds_out = xr.Dataset(data_vars=None, coords=out.coords, attrs=ds.attrs)\n",
    "ds_out[out.name] = out\n",
    "\n",
    "output_file = output_folder / \"test_tx_max.nc\"\n",
    "ds_out.to_netcdf(output_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "keep_output": false
   },
   "source": [
    "#### Loading the data in memory\n",
    "If the dataset is relatively small, it might be more efficient to simply load the data into the memory and use  numpy arrays instead of dask arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "ds4 = ds3.load()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "keep_output": false
   },
   "source": [
    "## Unit handling in `xclim`\n",
    "\n",
    "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",
    "The [Units Handling](units.ipynb) page goes more into detail on how unit conversion can easily be done.\n",
    "\n",
    "In the example below, we compute weekly total precipitation in mm using inputs of `mm/s` and `mm/d`. As we can see, the output is identical."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "# Compute with the original mm s-1 data\n",
    "out1 = xclim.atmos.precip_accumulation(ds4.pr, freq=\"MS\")\n",
    "# Create a copy of the data converted to mm d-1\n",
    "pr_mmd = ds4.pr * 3600 * 24\n",
    "pr_mmd.attrs[\"units\"] = \"mm d-1\"\n",
    "out2 = xclim.atmos.precip_accumulation(pr_mmd, freq=\"MS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "out1.plot(label=\"From mm s-1\", linestyle=\"-\")\n",
    "out2.plot(label=\"From mm d-1\", linestyle=\"none\", marker=\"o\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "keep_output": false
   },
   "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. Alternatively, users can even provide a threshold in Kelvin (``\"293.15 K\"``, if they really wanted to)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": false
   },
   "outputs": [],
   "source": [
    "# Create a copy of the data converted to C\n",
    "tasmax_C = ds4.tasmax - 273.15\n",
    "tasmax_C.attrs[\"units\"] = \"C\"\n",
    "\n",
    "# Using Kelvin data, threshold in Celsius\n",
    "out1 = xclim.atmos.tx_days_above(ds4.tasmax, thresh=\"20 C\", freq=\"MS\")\n",
    "\n",
    "# Using Celsius data\n",
    "out2 = xclim.atmos.tx_days_above(tasmax_C, thresh=\"20 C\", freq=\"MS\")\n",
    "\n",
    "# Using Celsius but with threshold in Kelvin\n",
    "out3 = xclim.atmos.tx_days_above(tasmax_C, 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=\"degC and degC\", marker=\"s\", markersize=10, linestyle=\"none\")\n",
    "out3.plot(label=\"degC and K\", marker=\"o\", linestyle=\"none\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Spatially varying thresholds\n",
    "Thresholds can also be passed as DataArrays instead of single scalar values, allowing the computation to depend on one or more non-temporal dimensions. The `units` attribute must be set.\n",
    "\n",
    "Going back to the initial `ds`, we'll subset it and compute the length of the heat wave according to thresholds that vary along the latitude and longitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds5 = ds.sel(time=slice(\"1950\", \"1960\"), lat=slice(50, 46), lon=slice(-75, -71))\n",
    "ds5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The tasmin threshold is 7°C for the northern half of the domain and 11°C for the southern half.\n",
    "# (notice that the lat coordinate is in decreasing order : from north to south)\n",
    "thresh_tasmin = xr.DataArray(\n",
    "    [7] * 24 + [11] * 24, dims=(\"lat\",), coords={\"lat\": ds5.lat}, attrs={\"units\": \"°C\"}\n",
    ")\n",
    "# The tasmax threshold is 17°C for the western half of the domain and 21°C for the eastern half.\n",
    "thresh_tasmax = xr.DataArray(\n",
    "    [17] * 24 + [21] * 24, dims=(\"lon\",), coords={\"lon\": ds5.lon}, attrs={\"units\": \"°C\"}\n",
    ")\n",
    "\n",
    "out_hw2d = xclim.atmos.heat_wave_total_length(\n",
    "    tasmin=ds5.tasmin,\n",
    "    tasmax=ds5.tasmax,\n",
    "    thresh_tasmin=thresh_tasmin,\n",
    "    thresh_tasmax=thresh_tasmax,\n",
    "    freq=\"YS\",\n",
    "    window=3,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final map for year 1958, shows clear jumps across the 4 quadrants, which was expected with our space-dependent thresholds. Notice also how the `long_name` (printed on the colorbar label) mentions that the threshold comes from \"an array\". This imprecise metadata is a consequence of using `DataArray`-derived thresholds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_hw2d.sel(time=\"1958\").plot()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}