Basic Usage¶

Climate indicator computations¶

xclim is a library of climate indicators that operate on xarray objects. Indicators perform health checks on input data, converts units as needed, assign nans when input data is missing, and format outputs according to the Climate and Forecast (CF) convention. As the list of indicators has grown quite large, indicators are accessed through their realm (xclim.indicators.atmos, xclim.indicators.land and xclim.indicators.seaIce, plus the special xclim.indicators.convert and xclim.indicators.generic) to help browsing indicators by the domain they apply to. Some more indicators are in project/institution-specific collections : xclim.indicators.cf holds the indicators defined by Clix-Meta, xclim.indicators.icclim for ECAD/ICCLIM and xclim.indicators.anuclim for Anuclim. See Climate Indicators for a list of available indicators.

Indicators should not be confused with compute functions, which define the algorithmic layer of each indicator. Those functions perform no checks beyond units compliance, and should be considered as low-level functions. See the respective documentation on indicators and compute functions for more information. These “compute functions” were previously named “indices”.

[ ]:

from __future__ import annotations

import xarray as xr

import xclim
from xclim.testing import open_dataset

Indicators are calculated by opening a NetCDF-like file, accessing the variable of interest, and calling the appropriate indicator function, which returns a new xarray.Dataset with the result.

For this example, we’ll first open a demonstration dataset storing surface air temperature and compute the number of growing degree days (the sum of degrees above a certain threshold) at the yearly frequency. The time frequency parameter, here YS, is specified using pandas offset aliases. Note that some offsets might not be supported for non-standard calendars (e.g. 360_day), see the xarray.cftime_range documentation for details.

[ ]:

# Normally, we would use xarray to open a dataset, e.g.:
# ds = xr.open_dataset("your_file.nc")

# For this example, let's use a test dataset from xclim:
ds = open_dataset("ERA5/daily_surface_cancities_1990-1993.nc")
ds.tas

Indicators are called as if they were functions.

[ ]:

gdd = xclim.atmos.growing_degree_days(tas=ds.tas, thresh="10.0 degC", freq="YS")
gdd

growing_degree_days makes most sense with daily input, but could theoretically accept other source frequencies.

This indicator has one argument thresh that represents a temperature threshold. When a single (scalar) value must be given, xclim accepts a string of the value with units. It will not accept a simple number as units can’t be inferred. Usually, xclim will also accept DataArrays with units for this kind of argument, as long as xarray is able to broadcast and align it with the input variable.

Finally, like almost all indicators, the function takes a freq argument to specify over what time period it is computed. These are called “Offset Aliases” by Pandas and are the same as the resampling string arguments (in xr.DataArray.resample). Valid arguments are detailed in pandas docs (note that aliases involving “business” notions are not supported by xarray and thus could raise issues in xclim).

Attributes: xclim added a global history attribute to the dataset, as well as other attributes to the output variable: standard_name, cell_methods, long_name , description, with the latter two formatted to include some information about the inputs. Of course, xclim updated the units, but also added units_metadata, required by the CF-Conventions to specify how temperature units should be interpreted. Finally, it added an xclim_description attribute that summarizes the computation done, similar to the history but without all the datetime, version, and provenance tracking.

Units handling paradigm¶

Indicators and compute functions are written in order to be flexible as to the sampling frequency and units of the data. You can use growing_degree_days on, for example, the 6-hourly data, but the output will then be in degree-hour units ("K h"). Moreover, all units, even when untouched by the calculation, will be reformatted into a CF-compliant symbol format. This behaviour was chosen to ensure consistency between all compute functions.

Indicators will always convert their output data into a specific unit, the one listed in [indicator].outputs, a list of metadata stores. For example here:

[ ]:

xclim.atmos.growing_degree_days.outputs[0]

This doesn’t apply to the compute functions. Very few compute functions will convert their output to specific units; Rather, it is the dimensionality that will be consistent on output.

Conventions¶

As you may have noticed, the growing_degree_days function above was not told along which dimension to operate. In xclim, the temporal dimension is always assumed to be named "time". All functions which reduce or compute over that dimension will expect that name. If you ever have another name in your data, you can simply rename it like:

ds = ds.rename(T="time")

For other names and attributes, xclim tries to follow different sets of conventions. In particular, input data should follow the CF conventions whenever possible for variable attributes. Variable names are usually the ones used in CMIP, when they exist. When a function needs to read the spatial coordinates, it will usually rely on CF conventions, through cf-xarray.

Indicator objects¶

Indicators are the core class of xclim. They act as compute functions with health checks and metadata attributes. Some checks involve:

Warning on missing input metadata
Identifying periods where missing data impacts the calculation and omit calculations for those periods. Those are called “missing methods” and are detailed in section Health checks.
Appending process history and maintaining the historical provenance of file metadata.
Writing Climate and Forecast Convention compliant metadata based on the variables and indicators calculated.

If we run the growing_degree_days indicator over a non-daily dataset, we’ll be warned that the input data is not daily. That is, running xclim.atmos.growing_degree_days(ds.air, thresh='10.0 degC', freq='MS') will fail with a ValidationError:

[ ]:

# Show that data is not at a daily time frequency

ds6h = xr.tutorial.load_dataset("air_temperature")
xr.infer_freq(ds6h.time)

[ ]:

gdd = xclim.atmos.growing_degree_days(tas=ds6h.air, thresh="10.0 degC", freq="MS")
gdd

Resampling to a daily frequency and running the same indicator succeeds, but we will still get warnings from the CF metadata checks.

[ ]:

daily_ds = ds6h.resample(time="D").mean(keep_attrs=True)
gdd = xclim.atmos.growing_degree_days(daily_ds.air, thresh="10.0 degC", freq="YS")
gdd

To suppress the CF validation warnings, we can set xclim to send these warnings to the log instead of raising a warning or an error. We also could set data_validation='warn' to be able to run the indicator on non-daily data. These options can be set globally or within a context manager with set_options.

The missing method which determines if a period should be considered missing or not can be controlled through the check_missing option, globally or contextually. The main missing methods also have options that can be modified.

[ ]:

with xclim.set_options(
    check_missing="pct",
    missing_options={"pct": dict(tolerance=0.1)},
    cf_compliance="log",
):
    # Change the missing method to "percent", instead of the default "any"
    # Set the tolerance to 10%, periods with more than 10% of missing data
    #     in the input will be masked in the output.
    gdd = xclim.atmos.growing_degree_days(daily_ds.air, thresh="10.0 degC", freq="MS")
gdd

Some indicators also expose time-selection arguments as **indexer keywords. This allows to compute on a subset of the time coordinates, for example only on a specific season, month, or between two dates in every year. It relies on the select_time function. Some indicators will simply select the time period and run the calculations, while others will smartly perform the selection at the right time, when the order of operation makes a difference. All will pass the indexer kwargs to the missing value handling, ensuring that the missing values outside the valid time period are not considered.

The next example computes the annual sum of growing degree days over 10 °C, but only considering days from the 1st of April to the 30th of September.

[ ]:

with xclim.set_options(cf_compliance="log"):
    gdd = xclim.atmos.growing_degree_days(tas=daily_ds.air, thresh="10 degC", freq="YS", date_bounds=("04-01", "09-30"))
gdd

xclim also allows us to call indicators using datasets and variable names. xarray’s keep_attrs option is used to decide if dataset (or node) attributes should be preserved and passed to the output.

[ ]:

with xclim.set_options(cf_compliance="log"):
    gdd = xclim.atmos.growing_degree_days(tas="air", thresh="10.0 degC", freq="MS", ds=daily_ds)

    # variable names default to xclim names, so we can even do this:
    renamed_daily_ds = daily_ds.rename(air="tas")
    gdd = xclim.atmos.growing_degree_days(thresh="10.0 degC", freq="MS", ds=renamed_daily_ds)
gdd

We can also get DataArrays as an output by disabling the as_dataset option.

[ ]:

with xclim.set_options(as_dataset=False, cf_compliance="log"):
    gdd_da = xclim.atmos.growing_degree_days(
        tas=daily_ds.air, thresh="10 degC", freq="YS", date_bounds=("04-01", "09-30")
    )
gdd_da

Finally, when passed a DataTree, xclim will map the computation over all nodes. It will skip empty nodes, but it requires that all non-empty nodes have all needed variables.

[ ]:

DT = xr.DataTree.from_dict({"/daily": daily_ds, "/daily_too": daily_ds, "/daily/again": daily_ds})
with xclim.set_options(cf_compliance="log"):
    gdd_dt = xclim.atmos.growing_degree_days(tas="air", thresh="10 degC", freq="YS", ds=DT)
gdd_dt

Graphics¶

xclim does not have specific functions to create graphics. However, it is built to ensure that indicators always have appropriate axis-related metadata that libraries like Matplotlib depend on to generate detailed and informative graphics.

This graphical functionality is entirely thanks to xarray, so the following examples are applicable to generic xarray.DataArray objects. For more examples, see the directions suggested by xarray’s plotting documentation

The xarray plot functions creates a histogram when the DataArray has 3 or more dimensions. In previous steps, xclim automatically filled the long_name and units attributes, which xarray uses to label the x-axis.

[ ]:

import matplotlib.pyplot as plt

print("long_name:", gdd.growing_degree_days.attrs["long_name"])
print("units:", gdd.growing_degree_days.attrs["units"])

gdd.growing_degree_days.plot()
plt.suptitle("Summary Statistics Histogram")
plt.show()

When the DataArray only has a time dimension, xarray plots a timeseries. In this case, xarray uses the long_name and units attributes provided by xclim to label the y-axis.

[ ]:

gdd.growing_degree_days.isel(lon=20, lat=10).plot()
plt.suptitle("Time Series at a Given Geographical Coordinate")
plt.show()

When the DataArray only has 2 dimensions, xarray plots a heatmap. In this case, xarray uses the long_name and units attributes provided by xclim to label the colorbar.

[ ]:

gdd.growing_degree_days.sel(time="2013-07-01").plot()
plt.suptitle("Spatial Pattern at a Specific Time Period")
plt.show()

Writing DataArrays and Datasets to disk¶

To save the data as a new NetCDF, use to_netcdf:

[ ]:

gdd.to_netcdf("monthly_growing_degree_days_data.nc")

It’s possible to save Dataset objects to other file formats. For more information see: xarray’s documentation