{
"cells": [
{
"cell_type": "markdown",
"id": "0",
"metadata": {},
"source": [
"# Climate change analysis of hydrological data"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1",
"metadata": {},
"outputs": [],
"source": [
"# Imports\n",
"from pathlib import Path\n",
"\n",
"import hvplot.xarray # noqa\n",
"import numpy as np\n",
"import pooch\n",
"import xarray as xr\n",
"import xclim\n",
"\n",
"import xhydro as xh\n",
"from xhydro.testing.helpers import deveraux\n",
"\n",
"D = deveraux()\n",
"\n",
"# Future streamflow file (1 file - Hydrotel driven by BCC-CSM-1.1(m))\n",
"streamflow_file = D.fetch(\"cc_indicators/streamflow_BCC-CSM1.1-m_rcp45.nc\")\n",
"\n",
"# Reference mean annual streamflow (QMOYAN) for 6 calibrations of Hydrotel\n",
"reference_files = D.fetch(\"cc_indicators/reference.zip\", pooch.Unzip())\n",
"\n",
"# Future deltas of QMOYAN (63 simulations x 6 calibrations of Hydrotel)\n",
"deltas_files = D.fetch(\"cc_indicators/deltas.zip\", pooch.Unzip())"
]
},
{
"cell_type": "markdown",
"id": "2",
"metadata": {},
"source": [
"While there is a vast array of analyses that can be performed to assess the impacts of climate change on hydrology, this notebook covers some of the most common steps:\n",
"\n",
"- Computing a list of relevant indicators over climatological periods.\n",
"- Computing future differences to assess the changes.\n",
"- Computing ensemble statistics to evaluate future changes and variability.\n",
"\n",
"INFO\n",
"\n",
"Several functions from the `xscen` library have been integrated into `xhydro` to simplify access for users, such as those in `xhydro.indicators` and `xhydro.cc`. This notebook will cover the basics, but for further details on these functions, please refer to the following resources:\n",
"\n",
"- [compute_indicators](https://xscen.readthedocs.io/en/latest/notebooks/2_getting_started.html#Computing-indicators)\n",
"- [climatological_op](https://xscen.readthedocs.io/en/latest/notebooks/2_getting_started.html#Climatological-operations)\n",
"- [compute_deltas](https://xscen.readthedocs.io/en/latest/notebooks/2_getting_started.html#Computing-deltas)\n",
"- [ensemble_statistics](https://xscen.readthedocs.io/en/latest/notebooks/2_getting_started.html#Ensemble-statistics)\n",
"\n",
"
\n"
]
},
{
"cell_type": "markdown",
"id": "3",
"metadata": {},
"source": [
"## Computing hydrological indicators over a given time period\n",
"\n",
"Hydrological indicators can be categorized into two main types:\n",
"\n",
"- Frequential indicators: These indicators describe hydrological events that occur at recurring intervals. They include metrics like the maximum 20-year flow (`Qmax20`) or the minimum 2-year 7-day average flow in summer (`Q7min2_summer`). The methodology for computing these indicators is covered in the [Local Frequency Analysis](../local_frequency_analysis.ipynb) notebook.\n",
"- Non-frequential indicators: These indicators do not explicitly describe recurrence, but rather absolute values or trends in hydrological variables. They include metrics like average yearly flow.\n",
"\n",
"Since frequential indicators are already covered in another example, this notebook will focus on the methodology for computing non-frequential indicators using `xhydro.indicators.compute_indicators`. This function is built on top of `xclim` and supports both predefined indicators, such as `xclim.indicator.land.doy_qmax`, as well as custom indicators created using `xclim.core.indicator.Indicator.from_dict`. The latter option can be quite complex—see the box below for more information. For advanced users, indicator construction can also be defined through a YAML file.\n",
"\n",
"The output of `xhydro.indicators.compute_indicators` is a dictionary, where each key represents the frequency of the requested indicators, following the `pandas` nomenclature. In our example, we will only use yearly data starting in January, so the frequency will be `YS-JAN`.\n",
"\n",
"\n",
" INFO\n",
"\n",
"Custom indicators in `xHydro` are built by following the YAML formatting required by `xclim`.\n",
"\n",
"A custom indicator built using `xclim.core.indicator.Indicator.from_dict` will need these elements:\n",
"\n",
"- \"data\": A dictionary with the following information:\n",
" - \"base\": The \"YAML ID\" obtained from [here](https://xclim.readthedocs.io/en/stable/indicators.html).\n",
" - \"input\": A dictionary linking the default xclim input to the name of your variable. Needed only if it is different. In the link above, they are the string following \"Uses:\".\n",
" - \"parameters\": A dictionary containing all other parameters for a given indicator. In the link above, the easiest way to access them is by clicking the link in the top-right corner of the box describing a given indicator.\n",
" - More entries can be used here, as described [in the xclim documentation](https://xclim.readthedocs.io/en/latest/api.html#yaml-file-structure) under \"identifier\".\n",
"- \"identifier\": A custom name for your indicator. This will be the name returned in the results.\n",
"- \"module\": Needed, but can be anything. To prevent an accidental overwriting of `xclim` indicators, it is best to use something different from: [\"atmos\", \"land\", \"generic\"].\n",
"\n",
"
\n"
]
},
{
"cell_type": "markdown",
"id": "4",
"metadata": {},
"source": [
"The example file used in this notebook is a daily time series of streamflow data, generated from the HYDROTEL hydrological model. This data is driven by bias-adjusted outputs from the BCC-CSM-1.1(m) climatological model (RCP4.5), spanning the years 1950 to 2100. For this example, the dataset includes data from just 2 stations. The function `xhydro.indicators.compute_indicators` can be used with any number of indicators. For this example, we will compute the mean annual flow and the mean summer-fall flow.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5",
"metadata": {},
"outputs": [],
"source": [
"ds = xr.open_dataset(streamflow_file).rename({\"streamflow\": \"q\"})\n",
"ds.q.hvplot(x=\"time\", grid=True, widget_location=\"bottom\", groupby=\"station\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6",
"metadata": {},
"outputs": [],
"source": [
"help(xh.indicators.compute_indicators)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7",
"metadata": {},
"outputs": [],
"source": [
"help(xclim.core.indicator.Indicator.from_dict)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8",
"metadata": {},
"outputs": [],
"source": [
"indicators = [\n",
" # 1st indicator: Mean annual flow\n",
" xclim.core.indicator.Indicator.from_dict(\n",
" data={\n",
" \"base\": \"stats\",\n",
" \"input\": {\"da\": \"q\"},\n",
" \"parameters\": {\"op\": \"mean\"},\n",
" },\n",
" identifier=\"QMOYAN\",\n",
" module=\"hydro\",\n",
" ),\n",
" # 2nd indicator: Mean summer-fall flow\n",
" xclim.core.indicator.Indicator.from_dict(\n",
" data={\n",
" \"base\": \"stats\",\n",
" \"input\": {\"da\": \"q\"},\n",
" \"parameters\": {\"op\": \"mean\", \"indexer\": {\"month\": [6, 7, 8, 9, 10, 11]}},\n",
" }, # The indexer is used to restrict available data to the relevant months only\n",
" identifier=\"QMOYEA\",\n",
" module=\"hydro\",\n",
" ),\n",
"]\n",
"\n",
"# Call compute_indicators\n",
"dict_indicators = xh.indicators.compute_indicators(ds, indicators=indicators)\n",
"\n",
"dict_indicators"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9",
"metadata": {},
"outputs": [],
"source": [
"dict_indicators[\"YS-JAN\"].QMOYAN.hvplot(\n",
" x=\"time\", grid=True, widget_location=\"bottom\", groupby=\"station\"\n",
")"
]
},
{
"cell_type": "markdown",
"id": "10",
"metadata": {},
"source": [
"The next step is to compute averages over climatological periods. This can be done using the `xhydro.cc.climatological_op` function.\n",
"\n",
"If the indicators themselves are not relevant to your analysis and you only need the climatological averages, you can directly use `xhydro.cc.produce_horizon` instead of combining `xhydro.indicators.compute_indicators` with `xhydro.cc.climatological_op`. The key advantage of `xhydro.cc.produce_horizon` is that it eliminates the `time` axis, replacing it with a `horizon` dimension that represents a slice of time. This is particularly useful when computing indicators with different output frequencies. An example of this approach is provided in the [Use Case Example](../use_case.ipynb).\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "11",
"metadata": {},
"outputs": [],
"source": [
"help(xh.cc.climatological_op)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "12",
"metadata": {},
"outputs": [],
"source": [
"# Call climatological_op. Here we don't need 'time' anymore, so we can use horizons_as_dim=True\n",
"ds_avg = xh.cc.climatological_op(\n",
" dict_indicators[\"YS-JAN\"],\n",
" op=\"mean\",\n",
" periods=[[1981, 2010], [2011, 2040], [2041, 2070], [2071, 2100]],\n",
" min_periods=29,\n",
" horizons_as_dim=True,\n",
" rename_variables=False,\n",
").drop_vars([\"time\"])\n",
"ds_avg"
]
},
{
"cell_type": "markdown",
"id": "13",
"metadata": {},
"source": [
"Once the averages over time periods have been computed, calculating the differences between future and past values is straightforward. Simply call `xhydro.cc.compute_deltas` to perform this calculation.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "14",
"metadata": {},
"outputs": [],
"source": [
"help(xh.cc.compute_deltas)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "15",
"metadata": {},
"outputs": [],
"source": [
"ds_deltas = xh.cc.compute_deltas(\n",
" ds_avg, reference_horizon=\"1981-2010\", kind=\"%\", rename_variables=False\n",
")\n",
"ds_deltas"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "16",
"metadata": {},
"outputs": [],
"source": [
"# Show the results as Dataframes\n",
"print(\"30-year averages\")\n",
"display(ds_avg.QMOYAN.isel(station=0).to_dataframe())\n",
"print(\"Deltas\")\n",
"display(ds_deltas.QMOYAN.isel(station=0).to_dataframe())"
]
},
{
"cell_type": "markdown",
"id": "17",
"metadata": {},
"source": [
"## Ensemble statistics\n",
"\n",
"In a real-world application, the steps outlined so far would need to be repeated for all available hydroclimatological simulations. For this example, we will work with a subset of pre-computed deltas from the RCP4.5 simulations used in the 2022 Hydroclimatic Atlas of Southern Quebec.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "18",
"metadata": {},
"outputs": [],
"source": [
"ds_dict_deltas = {}\n",
"for f in deltas_files:\n",
" id = Path(f).stem\n",
" ds_dict_deltas[id] = xr.open_dataset(f)\n",
"\n",
"print(f\"Loaded data from {len(ds_dict_deltas)} simulations\")"
]
},
{
"cell_type": "markdown",
"id": "19",
"metadata": {},
"source": [
"It is considered good practice to use multiple climate models when performing climate change analyses, especially since the impacts on the hydrological cycle can be nonlinear. Once multiple hydrological simulations are completed and ready for analysis, you can use `xhydro.cc.ensemble_stats` to access a variety of functions available in `xclim.ensemble`, such as calculating ensemble quantiles or assessing the agreement on the sign of change.\n",
"\n",
"### Weighting simulations\n",
"\n",
"When the ensemble of climate models is heterogeneous—such as when one model provides more simulations than others—it is recommended to weight the results accordingly. While this functionality is not currently available directly through `xhydro` (as it expects metadata specific to `xscen` workflows), the `xscen.generate_weights` function can help create an approximation of the weights based on available metadata.\n",
"\n",
"The following attributes are required for the function to work properly:\n",
"\n",
"- `'cat:source'` in all datasets\n",
"- `'cat:driving_model'` in regional climate models\n",
"- `'cat:institution'` in all datasets (if `independence_level='institution'`)\n",
"- `'cat:experiment'` in all datasets (if `split_experiments=True`)\n",
"\n",
"The `xscen.generate_weights` function offers three possible independence levels:\n",
"\n",
"- `model` (1 Model - 1 Vote): This assigns a total weight of 1 to all unique combinations of `'cat:source'` and `'cat:driving_model'`.\n",
"- `GCM` (1 GCM - 1 Vote): This assigns a total weight of 1 to all unique global climate models (GCMs), effectively averaging together regional climate simulations that originate from the same driving model.\n",
"- `institution` (1 institution - 1 Vote): This assigns a total weight of 1 to all unique `'cat:institution'` values.\n",
"\n",
"In all cases, the \"total weight of 1\" is not distributed equally between the involved simulations. The function will attempt to respect the model genealogy when distributing the weights. For example, if an institution has produced 4 simulations from Model A and 1 simulation from Model B, using `independence_level='institution'` would result in a weight of 0.125 for each Model A run and 0.5 for the single Model B run.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "20",
"metadata": {},
"outputs": [],
"source": [
"import xscen\n",
"\n",
"independence_level = \"model\" # 1 Model - 1 Vote\n",
"weights = xscen.generate_weights(ds_dict_deltas, independence_level=\"model\")\n",
"\n",
"# Show the results. We multiply by 6 for the showcase here simply because there are 6 hydrological platforms in the results.\n",
"weights.where(weights.realization.str.contains(\"LN24HA\"), drop=True) * 6"
]
},
{
"cell_type": "markdown",
"id": "21",
"metadata": {},
"source": [
"### Ensemble statistics with deterministic reference data\n",
"\n",
"In most cases, you will have deterministic data for the reference period. This means that, for a given location, the 30-year average for a specific indicator is represented by a single value.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "22",
"metadata": {},
"outputs": [],
"source": [
"# The Hydrological Portrait produces probabilistic estimates, but we'll take the 50th percentile to fake deterministic data\n",
"ref = xr.open_dataset(reference_files[0]).sel(percentile=50).drop_vars(\"percentile\")"
]
},
{
"cell_type": "markdown",
"id": "23",
"metadata": {},
"source": [
"Given that biases may still persist in climate simulations even after bias adjustment, which can impact hydrological modeling, we need to employ a perturbation technique to combine data over the reference period with climate simulations. This is particularly important in hydrology, where nonlinear interactions between climate and hydrological indicators can be significant. Multiple other methodologies exist for combining observed and simulated data, but comparing various approaches goes beyond the scope of this example.\n",
"\n",
"The perturbation technique involves calculating ensemble percentiles on the deltas and then applying those percentiles to the reference dataset. For this example, we'll compute the 10th, 25th, 50th, 75th, and 90th percentiles of the ensemble, as well as the agreement on the sign of the change, using the `xhydro.cc.ensemble_stats` function.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "24",
"metadata": {},
"outputs": [],
"source": [
"help(xh.cc.ensemble_stats)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "25",
"metadata": {},
"outputs": [],
"source": [
"statistics = {\n",
" \"ensemble_percentiles\": {\"values\": [10, 25, 50, 75, 90], \"split\": False},\n",
" \"robustness_fractions\": {\"test\": None},\n",
"}\n",
"\n",
"ens_stats = xh.cc.ensemble_stats(ds_dict_deltas, statistics, weights=weights)\n",
"ens_stats"
]
},
{
"cell_type": "markdown",
"id": "26",
"metadata": {},
"source": [
"This results in a large amount of data with many unique variables. To simplify the results, we'll compute three new statistics:\n",
"\n",
"- The median change.\n",
"- The interquartile range of the change.\n",
"- The agreement between models using the IPCC categories.\n",
"\n",
"The last statistic is slightly more complex. For more details on the categories of agreement for the sign of change, refer to the technical summary in \"Climate Change 2021 – The Physical Science Basis: Working Group I Contribution to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change\", [Cross-Chapter Box 1](https://www.cambridge.org/core/books/climate-change-2021-the-physical-science-basis/atlas/24E1C016DBBE4725BDFBC343695DE7DB). \n",
"\n",
"To compute this, you can use the results produced by `robustness_fractions`, but it needs a call to the function `xclim.ensembles.robustness_categories`. The thresholds and operations require two entries: the first is related to the significance test, and the second refers to the percentage of simulations showing a positive delta. For example, \"Agreement towards increase\" is met if more than 66% of simulations show a significant change, and 80% of simulations see a positive change.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "27",
"metadata": {},
"outputs": [],
"source": [
"out = xr.Dataset()\n",
"\n",
"out[\"QMOYAN_median\"] = ens_stats[\"QMOYAN\"].sel(percentiles=50)\n",
"out[\"QMOYAN_iqr\"] = ens_stats[\"QMOYAN\"].sel(percentiles=75) - ens_stats[\"QMOYAN\"].sel(\n",
" percentiles=25\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "28",
"metadata": {},
"outputs": [],
"source": [
"from xclim.ensembles import robustness_categories\n",
"\n",
"categories = [\n",
" \"Agreement towards increase\",\n",
" \"Agreement towards decrease\",\n",
" \"Conflicting signals\",\n",
" \"No change or robust signal\",\n",
"]\n",
"thresholds = [[0.66, 0.8], [0.66, 0.2], [0.66, 0.8], [0.66, np.nan]]\n",
"ops = [[\">=\", \">=\"], [\">=\", \"<=\"], [\">=\", \"<\"], [\"<\", None]]\n",
"\n",
"out[\"QMOYAN_robustness_categories\"] = robustness_categories(\n",
" changed_or_fractions=ens_stats[\"QMOYAN_changed\"],\n",
" agree=ens_stats[\"QMOYAN_positive\"],\n",
" categories=categories,\n",
" thresholds=thresholds,\n",
" ops=ops,\n",
")"
]
},
{
"cell_type": "markdown",
"id": "29",
"metadata": {},
"source": [
"Finally, using a perturbation method, future values for QMOYAN can be obtained by multiplying the reference indicator with the percentiles of the ensemble deltas."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "30",
"metadata": {},
"outputs": [],
"source": [
"out[\"QMOYAN_projected\"] = ref.QMOYAN * (1 + ens_stats.QMOYAN / 100)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "31",
"metadata": {},
"outputs": [],
"source": [
"out"
]
},
{
"cell_type": "markdown",
"id": "32",
"metadata": {},
"source": [
"### Ensemble statistics with probabilistic reference data\n",
"\n",
"This method is similar to the previous section, but it applies to cases like the [Hydrological Atlas of Southern Quebec](https://cehq.gouv.qc.ca/atlas-hydroclimatique/) or results from the [Optimal Interpolation](../optimal_interpolation.ipynb) notebook, where hydrological indicators for the historical period are represented by a probability density function (PDF) rather than a single discrete value. In such cases, the ensemble percentiles cannot simply be multiplied by the reference value.\n",
"\n",
"In this example, instead of a single value, `QMOYAN` is represented by 21 percentiles that capture the uncertainty surrounding this statistic. Similar to the future simulations, we also have 6 hydrological platforms to consider.\n",
"\n",
" WARNING\n",
"\n",
"In these cases, the percentiles in `ref` represent uncertainty (e.g., related to hydrological modeling or input data uncertainty), not interannual variability. At this stage, the temporal average should already have been calculated.\n",
"\n",
"
\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "33",
"metadata": {},
"outputs": [],
"source": [
"ref = xr.open_mfdataset(reference_files, combine=\"nested\", concat_dim=\"platform\")\n",
"\n",
"ref"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "34",
"metadata": {},
"outputs": [],
"source": [
"# This can also be represented as a cumulative distribution function (CDF)\n",
"import matplotlib.pyplot as plt\n",
"\n",
"for platform in ref.platform:\n",
" plt.plot(\n",
" ref.QMOYAN.isel(station=0).sel(platform=platform),\n",
" ref.QMOYAN.percentile / 100,\n",
" \"grey\",\n",
" )\n",
" plt.xlabel(\"Mean annual flow (m³/s)\")\n",
" plt.ylabel(\"Probability\")\n",
" plt.title(\"CDF for QMOYAN @ ABIT00057 \\nEach line is an hydrological platform\")"
]
},
{
"cell_type": "markdown",
"id": "35",
"metadata": {},
"source": [
"Due to their probabilistic nature, the historical reference values cannot be easily combined with the future deltas. To address this, the `xhydro.cc.weighted_random_sampling` and `xhydro.cc.sampled_indicators` functions have been designed. Together, these functions will:\n",
"\n",
"1. Sample 'n' values from the historical distribution, in accordance with the 'percentile' dimension.\n",
"2. Sample 'n' values from the delta distribution, using the provided weights.\n",
"3. Create the future distribution by applying the sampled deltas to the sampled historical distribution element-wise.\n",
"4. Compute the percentiles of the future distribution.\n",
"\n",
"First, we will sample within the reference dataset to combine the results of the 6 hydrological platforms together."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "36",
"metadata": {},
"outputs": [],
"source": [
"help(xh.cc.weighted_random_sampling)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "37",
"metadata": {},
"outputs": [],
"source": [
"deltas = xclim.ensembles.create_ensemble(ds_dict_deltas)\n",
"\n",
"hist_dist = xh.cc.weighted_random_sampling(\n",
" ds=ref,\n",
" include_dims=[\"platform\"],\n",
" n=10000,\n",
" seed=0,\n",
")\n",
"\n",
"hist_dist"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "38",
"metadata": {},
"outputs": [],
"source": [
"# Let's show how the historical distribution was sampled and reconstructed\n",
"def _make_cdf(ds, bins):\n",
" count, bins_count = np.histogram(ds.QMOYAN.isel(station=0), bins=bins)\n",
" pdf = count / sum(count)\n",
" return bins_count, np.cumsum(pdf)\n",
"\n",
"\n",
"# Barplot\n",
"plt.subplot(2, 1, 1)\n",
"uniquen = np.unique(hist_dist.QMOYAN.isel(station=0), return_counts=True)\n",
"plt.bar(uniquen[0], uniquen[1], width=0.01, color=\"k\")\n",
"plt.ylabel(\"Number of instances\")\n",
"plt.title(\"Sampling within the historical distribution\")\n",
"\n",
"# CDF\n",
"plt.subplot(2, 1, 2)\n",
"for i, platform in enumerate(ref.platform):\n",
" plt.plot(\n",
" ref.QMOYAN.isel(station=0).sel(platform=platform),\n",
" ref.percentile / 100,\n",
" \"grey\",\n",
" label=\"CDFs from the percentiles\" if i == 0 else None,\n",
" )\n",
"bc, c = _make_cdf(hist_dist, bins=50)\n",
"plt.plot(bc[1:], c, \"r\", label=f\"Sampled historical CDF (n={10000})\", linewidth=3)\n",
"plt.ylabel(\"Probability\")\n",
"plt.xlabel(\"QMOYAN (m³/s)\")\n",
"plt.legend()\n",
"\n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"id": "39",
"metadata": {},
"source": [
"We can do the same for the deltas. Since `weights` already contains all dimensions that we want to sample from, we don't need `include_dims` here."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "40",
"metadata": {},
"outputs": [],
"source": [
"delta_dist = xh.cc.weighted_random_sampling(\n",
" ds=deltas,\n",
" weights=weights,\n",
" n=10000,\n",
" seed=0,\n",
")\n",
"\n",
"delta_dist"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "41",
"metadata": {},
"outputs": [],
"source": [
"# Then, let's show how the deltas were sampled, for the last horizon\n",
"plt.subplot(2, 1, 1)\n",
"uniquen = np.unique(delta_dist.QMOYAN.isel(station=0, horizon=-1), return_counts=True)\n",
"plt.bar(uniquen[0], uniquen[1], width=0.25, color=\"k\")\n",
"plt.ylabel(\"Number of instances\")\n",
"plt.title(\"Sampling within the historical distribution\")\n",
"\n",
"plt.subplot(2, 1, 2)\n",
"bc, c = _make_cdf(delta_dist, bins=100)\n",
"plt.plot(bc[1:], c, \"k\", label=f\"Sampled deltas CDF (n={10000})\", linewidth=3)\n",
"plt.ylabel(\"Probability\")\n",
"plt.xlabel(\"Deltas (%)\")\n",
"plt.legend()\n",
"\n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"id": "42",
"metadata": {},
"source": [
"Once the two distributions have been acquired, `xhydro.cc.sampled_indicators` can be used to combine them element-wise and reconstruct a future distribution. The resulting distribution will possess the unique dimensions from both datasets. Here, this means that we get a reconstructed distribution for each future horizon."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "43",
"metadata": {},
"outputs": [],
"source": [
"help(xh.cc.sampled_indicators)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "44",
"metadata": {},
"outputs": [],
"source": [
"fut_dist, fut_pct = xh.cc.sampled_indicators(\n",
" ds_dist=hist_dist,\n",
" deltas_dist=delta_dist,\n",
" delta_kind=\"percentage\",\n",
" percentiles=ref.percentile,\n",
")\n",
"\n",
"fut_dist"
]
},
{
"cell_type": "markdown",
"id": "45",
"metadata": {},
"source": [
"Since we used the `percentiles` argument, it also computed a series of percentiles."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "46",
"metadata": {},
"outputs": [],
"source": [
"fut_pct"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "47",
"metadata": {},
"outputs": [],
"source": [
"# The distributions themselves can be used to create boxplots and compare the historical distribution to the future ones.\n",
"plt.boxplot(\n",
" [\n",
" hist_dist.QMOYAN.isel(station=0),\n",
" fut_dist.QMOYAN.isel(station=0, horizon=0),\n",
" fut_dist.QMOYAN.isel(station=0, horizon=1),\n",
" fut_dist.QMOYAN.isel(station=0, horizon=2),\n",
" ],\n",
" labels=[\"Historical\", \"2011-2040\", \"2041-2070\", \"2071-2100\"],\n",
")\n",
"\n",
"plt.ylabel(\"Mean summer flow (m³/s)\")\n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"id": "48",
"metadata": {},
"source": [
"The same statistics as before can also be computed by using the 10,000 samples within `delta_dist`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "49",
"metadata": {},
"outputs": [],
"source": [
"# The same statistics as before can also be computed by using delta_dist\n",
"delta_dist = delta_dist.rename({\"sample\": \"realization\"}) # xclim compatibility\n",
"ens_stats = xh.cc.ensemble_stats(delta_dist, statistics)\n",
"\n",
"out_prob = xr.Dataset()\n",
"out_prob[\"QMOYAN_median\"] = ens_stats[\"QMOYAN\"].sel(percentiles=50)\n",
"out_prob[\"QMOYAN_iqr\"] = ens_stats[\"QMOYAN\"].sel(percentiles=75) - ens_stats[\n",
" \"QMOYAN\"\n",
"].sel(percentiles=25)\n",
"out_prob[\"QMOYAN_robustness_categories\"] = robustness_categories(\n",
" changed_or_fractions=ens_stats[\"QMOYAN_changed\"],\n",
" agree=ens_stats[\"QMOYAN_positive\"],\n",
" categories=categories,\n",
" thresholds=thresholds,\n",
" ops=ops,\n",
")\n",
"\n",
"out_prob"
]
}
],
"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.11.14"
}
},
"nbformat": 4,
"nbformat_minor": 5
}