9. Analyse de l’impact des changements climatiques sur l’hydrologie¶
[1]:
# Imports
from pathlib import Path
import hvplot.xarray # noqa
import numpy as np
import pooch
import xarray as xr
import xclim
import xhydro as xh
from xhydro.testing.helpers import deveraux
D = deveraux()
# Future streamflow file (1 file - Hydrotel driven by BCC-CSM-1.1(m))
streamflow_file = D.fetch("cc_indicators/streamflow_BCC-CSM1.1-m_rcp45.nc")
# Reference mean annual streamflow (QMOYAN) for 6 calibrations of Hydrotel
reference_files = D.fetch("cc_indicators/reference.zip", pooch.Unzip())
# Future deltas of QMOYAN (63 simulations x 6 calibrations of Hydrotel)
deltas_files = D.fetch("cc_indicators/deltas.zip", pooch.Unzip())
ERROR 1: PROJ: proj_create_from_database: Open of /home/docs/checkouts/readthedocs.org/user_builds/xhydro-fr/conda/latest/share/proj failed
/home/docs/checkouts/readthedocs.org/user_builds/xhydro-fr/conda/latest/lib/python3.12/site-packages/xclim/sdba/__init__.py:22: FutureWarning: The SDBA submodule is in the process of being split from `xclim` in order to facilitate development and effective maintenance of the SDBA utilities. The `xclim.sdba` functionality will change in the future. For more information, please visit https://xsdba.readthedocs.io/en/latest/.
warnings.warn(
Bien qu’il existe une vaste gamme d’analyses qui peuvent être effectuées pour évaluer les impacts du changement climatique sur l’hydrologie, ce Notebook couvre quelques-unes des étapes les plus courantes :
Calculer une liste d’indicateurs pertinents sur des périodes climatologiques.
Calculer les différences futures pour évaluer les changements.
Calculer les statistiques d’ensemble pour évaluer les changements et la variabilité futurs.
INFO
Plusieurs fonctions de la librairie xscen ont été intégrées dans xhydro pour en simplifier l’accès aux utilisateurs, comme celles dans xhydro.indicators et xhydro.cc. Ce Notebook couvrira les bases, mais pour plus de détails sur ces fonctions, veuillez consulter les ressources suivantes :
9.1. Calcul des indicateurs hydrologiques sur une période donnée¶
Les indicateurs hydrologiques peuvent être classés en deux types principaux :
Indicateurs fréquentiels : Ces indicateurs décrivent les événements hydrologiques qui se produisent à intervalles réguliers. Ils incluent des métriques comme le débit maximal sur 20 ans (
Qmax20) ou le débit minimum moyenné sur 7 jours de récurrence 2 ans en été (Q7min2_summer). La méthodologie de calcul de ces indicateurs est abordée dans le Notebook Analyses fréquencielles locales.Indicateurs non fréquentiels : Ces indicateurs ne décrivent pas explicitement la récurrence, mais plutôt des valeurs absolues ou des tendances dans les variables hydrologiques. Ils incluent des métriques comme le débit moyen annuel.
Étant donné que les indicateurs fréquentiels sont déjà couverts dans un autre exemple, ce Notebook se concentrera sur la méthodologie de calcul des indicateurs non fréquentiels en utilisant xhydro.indicators.compute_indicators. Cette fonction est construite sur la base de xclim et prend en charge à la fois des indicateurs prédéfinis, tels que xclim.indicator.land.doy_qmax, ainsi que des indicateurs personnalisés créés à l’aide de xclim.core.indicator.Indicator.from_dict. Cette dernière option peut être assez complexe—voir la boîte ci-dessous pour plus d’informations. Pour les utilisateurs avancés, la construction d’indicateurs peut également être définie via un fichier YAML.
La sortie de xhydro.indicators.compute_indicators est un dictionnaire, où chaque clé représente la fréquence des indicateurs demandés, suivant la nomenclature pandas. Dans notre exemple, nous n’utiliserons que les données annuelles commençant en janvier, donc la fréquence sera YS-JAN.
INFO
Les indicateurs personnalisés dans xHydro sont construits en suivant le formatage YAML requis par xclim.
Un indicateur personnalisé construit avec xclim.core.indicator.Indicator.from_dict nécessitera ces éléments :
« data » : Un dictionnaire avec les informations suivantes :
« base » : L“« ID YAML » obtenu depuis cette page.
« input » : Un dictionnaire reliant l’entrée xclim par défaut au nom de votre variable. Nécessaire uniquement si elle diffère. Dans le lien ci-dessus, ce sont les chaînes suivant « Uses: ».
« parameters » : Un dictionnaire contenant tous les autres paramètres pour un indicateur donné. Dans le lien ci-dessus, la manière la plus simple d’y accéder est de cliquer sur le lien dans le coin supérieur droit de la boîte décrivant un indicateur donné. »é.
Des entrées supplémentaires peuvent être utilisées ici, comme décrit dans la documentation xclim sous « identifier ».
« identifier » : Un nom personnalisé pour votre indicateur. Ce sera le nom retourné dans les résultats.
« module »: Nécessaire, mais peut être n’importe quoi. Pour éviter un écrasement accidentel des indicateurs
xclim, il est préférable d’utiliser quelque chose de différent de : [« atmos », « land », « generic »].
Le fichier d’exemple utilisé dans ce Notebook est une série temporelle quotidienne des débits, générée à partir du modèle hydrologique HYDROTEL. Ces données sont pilotées par des sorties post-traitées du modèle climatique BCC-CSM-1.1(m) (RCP4.5), couvrant les années 1950 à 2100. Pour cet exemple, le jeu de données comprend des données de seulement 2 stations. La fonction xhydro.indicators.compute_indicators peut être utilisée avec n’importe quel nombre d’indicateurs. Pour cet exemple, nous allons calculer le débit moyen annuel et le débit moyen été-automne.
[2]:
ds = xr.open_dataset(streamflow_file).rename({"streamflow": "q"})
ds.q.hvplot(x="time", grid=True, widget_location="bottom", groupby="station")
[2]:
[3]:
Help on function compute_indicators in module xscen.indicators:
compute_indicators(ds: xarray.core.dataset.Dataset, indicators: str | os.PathLike | collections.abc.Sequence[xclim.core.indicator.Indicator] | collections.abc.Sequence[tuple[str, xclim.core.indicator.Indicator]] | module, *, periods: list[str] | list[list[str]] | None = None, restrict_years: bool = True, to_level: str | None = 'indicators', rechunk_input: bool = False) -> dict
Calculate variables and indicators based on a YAML call to xclim.
The function cuts the output to be the same years as the inputs.
Hence, if an indicator creates a timestep outside the original year range (e.g. the first DJF for QS-DEC),
it will not appear in the output.
Parameters
----------
ds : xr.Dataset
Dataset to use for the indicators.
indicators : str | os.PathLike | Sequence[Indicator] | Sequence[tuple[str, Indicator]] | ModuleType
Path to a YAML file that instructs on how to calculate missing variables.
Can also be only the "stem", if translations and custom indices are implemented.
Can be the indicator module directly, or a sequence of indicators or a sequence of
tuples (indicator name, indicator) as returned by `iter_indicators()`.
periods : list of str or list of lists of str, optional
Either [start, end] or list of [start, end] of continuous periods over which to compute the indicators.
This is needed when the time axis of ds contains some jumps in time.
If None, the dataset will be considered continuous.
restrict_years : bool
If True, cut the time axis to be within the same years as the input.
This is mostly useful for frequencies that do not start in January, such as QS-DEC.
In that instance, `xclim` would start on previous_year-12-01 (DJF), with a NaN.
`restrict_years` will cut that first timestep. This should have no effect on YS and MS indicators.
to_level : str, optional
The processing level to assign to the output.
If None, the processing level of the inputs is preserved.
rechunk_input : bool
If True, the dataset is rechunked with :py:func:`flox.xarray.rechunk_for_blockwise`
according to the resampling frequency of the indicator. Each rechunking is done
once per frequency with :py:func:`xscen.utils.rechunk_for_resample`.
Returns
-------
dict
Dictionary (keys = timedeltas) with indicators separated by temporal resolution.
See Also
--------
xclim.indicators, xclim.core.indicator.build_indicator_module_from_yaml
[4]:
Help on method from_dict in module xclim.core.indicator:
from_dict(data: 'dict', identifier: 'str', module: 'str | None' = None) -> 'Indicator' class method of xclim.core.indicator.Indicator
Create an indicator subclass and instance from a dictionary of parameters.
Most parameters are passed directly as keyword arguments to the class constructor, except:
- "base" : A subclass of Indicator or a name of one listed in
:py:data:`xclim.core.indicator.registry` or
:py:data:`xclim.core.indicator.base_registry`. When passed, it acts as if
`from_dict` was called on that class instead.
- "compute" : A string function name translates to a
:py:mod:`xclim.indices.generic` or :py:mod:`xclim.indices` function.
Parameters
----------
data : dict
The exact structure of this dictionary is detailed in the submodule documentation.
identifier : str
The name of the subclass and internal indicator name.
module : str
The module name of the indicator. This is meant to be used only if the indicator
is part of a dynamically generated submodule, to override the module of the base class.
Returns
-------
Indicator
A new Indicator instance.
[5]:
indicators = [
# 1st indicator: Mean annual flow
xclim.core.indicator.Indicator.from_dict(
data={
"base": "stats",
"input": {"da": "q"},
"parameters": {"op": "mean"},
},
identifier="QMOYAN",
module="hydro",
),
# 2nd indicator: Mean summer-fall flow
xclim.core.indicator.Indicator.from_dict(
data={
"base": "stats",
"input": {"da": "q"},
"parameters": {"op": "mean", "indexer": {"month": [6, 7, 8, 9, 10, 11]}},
}, # The indexer is used to restrict available data to the relevant months only
identifier="QMOYEA",
module="hydro",
),
]
# Call compute_indicators
dict_indicators = xh.indicators.compute_indicators(ds, indicators=indicators)
dict_indicators
[5]:
{'YS-JAN': <xarray.Dataset> Size: 4kB
Dimensions: (station: 2, time: 151)
Coordinates:
lat (station) float32 8B ...
lon (station) float32 8B ...
drainage_area (station) float32 8B ...
* station (station) <U9 72B 'ABIT00057' 'ABIT00058'
* time (time) datetime64[ns] 1kB 1950-01-01 ... 2100-01-01
Data variables:
QMOYAN (station, time) float32 1kB nan nan nan ... 0.8027 0.743 nan
QMOYEA (station, time) float32 1kB nan nan nan ... 0.385 1.035
Attributes: (12/48)
institution: DEH (Direction de l'Expertise Hydrique)
institute_id: DEH
contact: atlas.hydroclimatique@environnement.gouv....
date_created: 2020-05-23
featureType: timeSeries
cdm_data_type: station
... ...
cat:bias_adjust_project: INFO-Crue-CM5
cat:mip_era: CMIP5
cat:activity: CMIP5
cat:experiment: rcp45
cat:member: r1i1p1
cat:bias_adjust_institution: OURANOS}
[6]:
dict_indicators["YS-JAN"].QMOYAN.hvplot(
x="time", grid=True, widget_location="bottom", groupby="station"
)
[6]:
La prochaine étape consiste à calculer des moyennes sur des périodes climatiques. Cela peut être fait en utilisant la fonction xhydro.cc.climatological_op.
Si les indicateurs eux-mêmes ne sont pas pertinents pour votre analyse et que vous n’avez besoin que des moyennes climatiques, vous pouvez directement utiliser xhydro.cc.produce_horizon au lieu de combiner xhydro.indicators.compute_indicators avec xhydro.cc.climatological_op. L’avantage clé de xhydro.cc.produce_horizon est qu’il élimine l’axe time, le remplaçant par une dimension horizon représentant une tranche de temps. Cela est particulièrement utile lors du calcul d’indicateurs avec des fréquences de sortie différentes. Un exemple de cette approche est fourni dans l”exemple d’application.
[7]:
help(xh.cc.climatological_op)
Help on function climatological_op in module xscen.aggregate:
climatological_op(ds: xarray.core.dataset.Dataset, *, op: str | dict = 'mean', window: int | None = None, min_periods: int | float | None = None, stride: int = 1, periods: list[str] | list[list[str]] | None = None, rename_variables: bool = True, to_level: str = 'climatology', horizons_as_dim: bool = False) -> xarray.core.dataset.Dataset
Perform an operation 'op' over time, for given time periods, respecting the temporal resolution of ds.
Parameters
----------
ds : xr.Dataset
Dataset to use for the computation.
op : str or dict
Operation to perform over time.
The operation can be any method name of xarray.core.rolling.DatasetRolling, 'linregress', or a dictionary.
If 'op' is a dictionary, the key is the operation name and the value is a dict of kwargs
accepted by the equivalent NumPy function. See the Notes for more information.
While other operations are technically possible, the following are recommended and tested:
['max', 'mean', 'median', 'min', 'std', 'sum', 'var', 'linregress'].
Operations beyond methods of xarray.core.rolling.DatasetRolling include:
- 'linregress' : Computes the linear regression over time, using
scipy.stats.linregress and employing years as regressors.
The output will have a new dimension 'linreg_param' with coordinates:
['slope', 'intercept', 'rvalue', 'pvalue', 'stderr', 'intercept_stderr'].
Only one operation per call is supported, so len(op)==1 if a dict.
window : int, optional
Number of years to use for the rolling operation.
If left at None and periods is given, window will be the size of the first period. Hence, if periods are of
different lengths, the shortest period should be passed first.
If left at None and periods is not given, the window will be the size of the input dataset.
min_periods : int or float, optional
For the rolling operation, minimum number of years required for a value to be computed.
If left at None and the xrfreq is either QS or AS and doesn't start in January,
min_periods will be one less than window.
Otherwise, if left at None, it will be deemed the same as 'window'.
If passed as a float value between 0 and 1, this will be interpreted as the floor of the fraction of the window size.
stride : int
Stride (in years) at which to provide an output from the rolling window operation.
periods : list of str or list of lists of str, optional
Either [start, end] or list of [start, end] of continuous periods to be considered.
This is needed when the time axis of ds contains some jumps in time.
If None, the dataset will be considered continuous.
rename_variables : bool
If True, '_clim_{op}' will be added to variable names.
to_level : str, optional
The processing level to assign to the output.
If None, the processing level of the inputs is preserved.
horizons_as_dim : bool
If True, the output will have 'horizon' and the frequency as 'month', 'season' or 'year' as
dimensions and coordinates. The 'time' coordinate will be unstacked to horizon and frequency dimensions.
Horizons originate from periods and/or windows and their stride in the rolling operation.
Returns
-------
xr.Dataset
Dataset with the results from the climatological operation.
Notes
-----
xarray.core.rolling.DatasetRolling functions do not support kwargs other than 'keep_attrs'. In order to pass
additional arguments to the operation, we instead use the 'reduce' method and pass the operation as a numpy function.
If possible, a function that handles NaN values will be used (e.g. op='mean' will use `np.nanmean`), as the
'min_periods' argument already decides how many NaN values are acceptable.
[8]:
# Call climatological_op. Here we don't need 'time' anymore, so we can use horizons_as_dim=True
ds_avg = xh.cc.climatological_op(
dict_indicators["YS-JAN"],
op="mean",
periods=[[1981, 2010], [2011, 2040], [2041, 2070], [2071, 2100]],
min_periods=29,
horizons_as_dim=True,
rename_variables=False,
).drop_vars(["time"])
ds_avg
[8]:
<xarray.Dataset> Size: 304B
Dimensions: (horizon: 4, station: 2)
Coordinates:
lat (station) float32 8B 49.53 49.53
lon (station) float32 8B -77.05 -77.02
drainage_area (station) float32 8B 60.76 54.61
* station (station) <U9 72B 'ABIT00057' 'ABIT00058'
* horizon (horizon) <U9 144B '1981-2010' '2011-2040' ... '2071-2100'
Data variables:
QMOYAN (horizon, station) float32 32B 0.8968 0.8041 ... 0.7745
QMOYEA (horizon, station) float32 32B 0.9335 0.8379 ... 0.6966
Attributes: (12/48)
institution: DEH (Direction de l'Expertise Hydrique)
institute_id: DEH
contact: atlas.hydroclimatique@environnement.gouv....
date_created: 2020-05-23
featureType: timeSeries
cdm_data_type: station
... ...
cat:bias_adjust_project: INFO-Crue-CM5
cat:mip_era: CMIP5
cat:activity: CMIP5
cat:experiment: rcp45
cat:member: r1i1p1
cat:bias_adjust_institution: OURANOSUne fois les moyennes sur les périodes de temps calculées, calculer les différences entre les valeurs futures et passées est simple. Il suffit d’appeler xhydro.cc.compute_deltas pour effectuer ce calcul.
[9]:
help(xh.cc.compute_deltas)
Help on function compute_deltas in module xscen.aggregate:
compute_deltas(ds: xarray.core.dataset.Dataset, reference_horizon: str | xarray.core.dataset.Dataset, *, kind: str | dict = '+', rename_variables: bool = True, to_level: str | None = 'deltas') -> xarray.core.dataset.Dataset
Compute deltas in comparison to a reference time period, respecting the temporal resolution of ds.
Parameters
----------
ds : xr.Dataset
Dataset to use for the computation.
reference_horizon : str or xr.Dataset
Either a YYYY-YYYY string corresponding to the 'horizon' coordinate of the reference period,
or a xr.Dataset containing the climatological mean.
kind : str or dict
['+', '/', '%'] Whether to provide absolute, relative, or percentage deltas.
Can also be a dictionary separated per variable name.
rename_variables : bool
If True, '_delta_YYYY-YYYY' will be added to variable names.
to_level : str, optional
The processing level to assign to the output.
If None, the processing level of the inputs is preserved.
Returns
-------
xr.Dataset
Returns a Dataset with the requested deltas.
[10]:
ds_deltas = xh.cc.compute_deltas(
ds_avg, reference_horizon="1981-2010", kind="%", rename_variables=False
)
ds_deltas
[10]:
<xarray.Dataset> Size: 304B
Dimensions: (station: 2, horizon: 4)
Coordinates:
lat (station) float32 8B 49.53 49.53
lon (station) float32 8B -77.05 -77.02
drainage_area (station) float32 8B 60.76 54.61
* station (station) <U9 72B 'ABIT00057' 'ABIT00058'
* horizon (horizon) <U9 144B '1981-2010' '2011-2040' ... '2071-2100'
Data variables:
QMOYAN (horizon, station) float32 32B 0.0 0.0 4.273 ... -3.64 -3.678
QMOYEA (horizon, station) float32 32B 0.0 0.0 ... -16.67 -16.86
Attributes: (12/48)
institution: DEH (Direction de l'Expertise Hydrique)
institute_id: DEH
contact: atlas.hydroclimatique@environnement.gouv....
date_created: 2020-05-23
featureType: timeSeries
cdm_data_type: station
... ...
cat:bias_adjust_project: INFO-Crue-CM5
cat:mip_era: CMIP5
cat:activity: CMIP5
cat:experiment: rcp45
cat:member: r1i1p1
cat:bias_adjust_institution: OURANOS[11]:
# Show the results as Dataframes
print("30-year averages")
display(ds_avg.QMOYAN.isel(station=0).to_dataframe())
print("Deltas")
display(ds_deltas.QMOYAN.isel(station=0).to_dataframe())
30-year averages
| lat | lon | drainage_area | station | QMOYAN | |
|---|---|---|---|---|---|
| horizon | |||||
| 1981-2010 | 49.529999 | -77.050003 | 60.759998 | ABIT00057 | 0.896769 |
| 2011-2040 | 49.529999 | -77.050003 | 60.759998 | ABIT00057 | 0.935091 |
| 2041-2070 | 49.529999 | -77.050003 | 60.759998 | ABIT00057 | 0.944270 |
| 2071-2100 | 49.529999 | -77.050003 | 60.759998 | ABIT00057 | 0.864129 |
Deltas
| lat | lon | drainage_area | station | QMOYAN | |
|---|---|---|---|---|---|
| horizon | |||||
| 1981-2010 | 49.529999 | -77.050003 | 60.759998 | ABIT00057 | 0.000000 |
| 2011-2040 | 49.529999 | -77.050003 | 60.759998 | ABIT00057 | 4.273284 |
| 2041-2070 | 49.529999 | -77.050003 | 60.759998 | ABIT00057 | 5.296827 |
| 2071-2100 | 49.529999 | -77.050003 | 60.759998 | ABIT00057 | -3.639797 |
9.2. Statistiques d’ensemble¶
Dans une application réelle, les étapes décrites jusqu’à présent devraient être répétées pour toutes les simulations hydroclimatiques disponibles. Pour cet exemple, nous travaillerons avec un sous-ensemble de deltas pré-calculés des simulations RCP4.5 utilisées dans l’édition 2022 de l’Atlas hydroclimatique du Québec méridional.
[12]:
ds_dict_deltas = {}
for f in deltas_files:
id = Path(f).stem
ds_dict_deltas[id] = xr.open_dataset(f)
print(f"Loaded data from {len(ds_dict_deltas)} simulations")
Loaded data from 378 simulations
Il est recommandé d’utiliser plusieurs modèles climatiques lors de l’analyse des changements climatiques, surtout que les impacts sur le cycle hydrologique peuvent être non linéaires. Une fois plusieurs simulations hydrologiques complétées et prêtes pour l’analyse, vous pouvez utiliser xhydro.cc.ensemble_stats pour accéder à une variété de fonctions disponibles dans xclim.ensemble, telles que le calcul des quantiles d’ensemble ou l’évaluation de l’accord sur le signe du changement.
9.2.1. Pondération des simulations¶
Lorsque l’ensemble de modèles climatiques est hétérogène—par exemple, lorsqu’un modèle fournit plus de simulations que les autres—il est recommandé de pondérer les résultats en conséquence. Bien que cette fonctionnalité ne soit pas actuellement disponible directement via xHydro (car elle attend des métadonnées spécifiques aux flux de travail xscen), la fonction xscen.generate_weights peut aider à créer une approximation des poids en fonction des métadonnées disponibles.
Les attributs suivants sont nécessaires pour que la fonction fonctionne correctement :
'cat:source'dans tous les jeux de données'cat:driving_model'dans les modèles climatiques régionaux'cat:institution'dans tous les jeux de données (siindependence_level='institution')'cat:experiment'dans tous les jeux de données (sisplit_experiments=True)
La fonction xscen.generate_weights propose trois niveaux d’indépendance possibles :
model(1 Modèle - 1 Vote) : Cela attribue un poids total de 1 à toutes les combinaisons uniques de'cat:source'et'cat:driving_model'.GCM(1 GCM - 1 Vote) : Cela attribue un poids total de 1 à tous les modèles climatiques globaux (GCM), en moyennant les simulations climatiques régionales provenant du même modèle pilote.institution(1 institution - 1 vote) : Cela attribue un poids total de 1 à toutes les valeurs uniques'cat:institution'.
Dans tous les cas, le « poids total de 1 » n’est pas réparti de manière égale entre les simulations impliquées. La fonction tentera de respecter la généalogie du modèle lors de la distribution des poids. Par exemple, si une institution a produit 4 simulations du Modèle A et 1 simulation du Modèle B, l’utilisation de independence_level='institution' donnerait un poids de 0,125 pour chaque simulation du Modèle A et 0,5 pour la simulation unique du Modèle B.
[13]:
import xscen
independence_level = "model" # 1 Model - 1 Vote
weights = xscen.generate_weights(ds_dict_deltas, independence_level="model")
# Show the results. We multiply by 6 for the showcase here simply because there are 6 hydrological platforms in the results.
weights.where(weights.realization.str.contains("LN24HA"), drop=True) * 6
[13]:
<xarray.DataArray (realization: 63)> Size: 504B
array([1. , 0.1 , 1. , 1. , 1. ,
1. , 0.5 , 0.33333333, 1. , 0.1 ,
1. , 0.33333333, 1. , 1. , 0.25 ,
1. , 1. , 1. , 0.33333333, 0.1 ,
1. , 1. , 1. , 0.1 , 0.2 ,
0.33333333, 0.5 , 0.25 , 1. , 0.1 ,
0.33333333, 0.1 , 0.5 , 0.5 , 0.25 ,
0.33333333, 1. , 0.33333333, 0.5 , 1. ,
0.1 , 0.1 , 1. , 0.2 , 1. ,
0.2 , 1. , 0.33333333, 0.33333333, 1. ,
0.5 , 0.1 , 0.33333333, 0.33333333, 1. ,
0.1 , 0.2 , 0.25 , 0.33333333, 0.2 ,
1. , 1. , 1. ])
Coordinates:
* realization (realization) <U112 28kB 'AtlasHydro2022_HYDROTEL_LN24HA_INF...9.2.2. Statistiques d’ensemble avec des données de référence déterministes¶
Dans la plupart des cas, vous disposerez de données déterministes pour la période de référence. Cela signifie que, pour un lieu donné, la moyenne sur 30 ans d’un indicateur spécifique est représentée par une seule valeur.
[14]:
# The Hydrological Portrait produces probabilistic estimates, but we'll take the 50th percentile to fake deterministic data
ref = xr.open_dataset(reference_files[0]).sel(percentile=50).drop_vars("percentile")
Étant donné que des biais peuvent encore persister dans les simulations climatiques même après ajustement des biais, ce qui peut affecter la modélisation hydrologique, nous devons employer une technique de perturbation pour combiner les données sur la période de référence avec les simulations climatiques. Cela est particulièrement important en hydrologie, où les interactions non linéaires entre le climat et les indicateurs hydrologiques peuvent être significatives. Plusieurs autres méthodologies existent pour combiner les données observées et simulées, mais la comparaison des différentes approches dépasse le cadre de cet exemple.
La technique de perturbation consiste à calculer les percentiles de l’ensemble sur les deltas, puis à appliquer ces percentiles au jeu de données de référence. Pour cet exemple, nous allons calculer les 10e, 25e, 50e, 75e et 90e percentiles de l’ensemble, ainsi que l’accord sur le signe du changement, en utilisant la fonction xhydro.cc.ensemble_stats.
[15]:
help(xh.cc.ensemble_stats)
Help on function ensemble_stats in module xscen.ensembles:
ensemble_stats(datasets: dict | list[str | os.PathLike] | list[xarray.core.dataset.Dataset] | list[xarray.core.dataarray.DataArray] | xarray.core.dataset.Dataset, statistics: dict, *, create_kwargs: dict | None = None, weights: xarray.core.dataarray.DataArray | None = None, common_attrs_only: bool = True, to_level: str = 'ensemble') -> xarray.core.dataset.Dataset
Create an ensemble and computes statistics on it.
Parameters
----------
datasets : dict or list of [str, os.PathLike, Dataset or DataArray], or Dataset
List of file paths or xarray Dataset/DataArray objects to include in the ensemble.
A dictionary can be passed instead of a list, in which case the keys are used as coordinates along the new
`realization` axis.
Tip: With a project catalog, you can do: `datasets = pcat.search(**search_dict).to_dataset_dict()`.
If a single Dataset is passed, it is assumed to already be an ensemble and will be used as is. The 'realization' dimension is required.
statistics : dict
xclim.ensembles statistics to be called. Dictionary in the format {function: arguments}.
If a function requires 'weights', you can leave it out of this dictionary and
it will be applied automatically if the 'weights' argument is provided.
See the Notes section for more details on robustness statistics, which are more complex in their usage.
create_kwargs : dict, optional
Dictionary of arguments for xclim.ensembles.create_ensemble.
weights : xr.DataArray, optional
Weights to apply along the 'realization' dimension. This array cannot contain missing values.
common_attrs_only : bool
If True, keeps only the global attributes that are the same for all datasets and generate new id.
If False, keeps global attrs of the first dataset (same behaviour as xclim.ensembles.create_ensemble)
to_level : str
The processing level to assign to the output.
Returns
-------
xr.Dataset
Dataset with ensemble statistics
Notes
-----
* The positive fraction in 'change_significance' and 'robustness_fractions' is calculated by
xclim using 'v > 0', which is not appropriate for relative deltas.
This function will attempt to detect relative deltas by using the 'delta_kind' attribute ('rel.', 'relative', '*', or '/')
and will apply 'v - 1' before calling the function.
* The 'robustness_categories' statistic requires the outputs of 'robustness_fractions'.
Thus, there are two ways to build the 'statistics' dictionary:
1. Having 'robustness_fractions' and 'robustness_categories' as separate entries in the dictionary.
In this case, all outputs will be returned.
2. Having 'robustness_fractions' as a nested dictionary under 'robustness_categories'.
In this case, only the robustness categories will be returned.
* A 'ref' DataArray can be passed to 'change_significance' and 'robustness_fractions', which will be used by xclim to compute deltas
and perform some significance tests. However, this supposes that both 'datasets' and 'ref' are still timeseries (e.g. annual means),
not climatologies where the 'time' dimension represents the period over which the climatology was computed. Thus,
using 'ref' is only accepted if 'robustness_fractions' (or 'robustness_categories') is the only statistic being computed.
* If you want to use compute a robustness statistic on a climatology, you should first compute the climatologies and deltas yourself,
then leave 'ref' as None and pass the deltas as the 'datasets' argument. This will be compatible with other statistics.
See Also
--------
xclim.ensembles._base.create_ensemble, xclim.ensembles._base.ensemble_percentiles,
xclim.ensembles._base.ensemble_mean_std_max_min,
xclim.ensembles._robustness.robustness_fractions, xclim.ensembles._robustness.robustness_categories,
xclim.ensembles._robustness.robustness_coefficient,
[16]:
statistics = {
"ensemble_percentiles": {"values": [10, 25, 50, 75, 90], "split": False},
"robustness_fractions": {"test": None},
}
ens_stats = xh.cc.ensemble_stats(ds_dict_deltas, statistics, weights=weights)
ens_stats
[16]:
<xarray.Dataset> Size: 820B
Dimensions: (station: 2, horizon: 3, percentiles: 5)
Coordinates:
lat (station) float32 8B 49.53 49.53
lon (station) float32 8B -77.05 -77.02
drainage_area (station) float32 8B 60.76 54.61
* station (station) <U9 72B 'ABIT00057' 'ABIT00058'
* horizon (horizon) <U9 108B '2021-2051' ... '2071-2100'
* percentiles (percentiles) int64 40B 10 25 50 75 90
Data variables:
QMOYAN (percentiles, station, horizon) float64 240B -3....
QMOYAN_changed (station, horizon) float64 48B 1.0 1.0 ... 1.0 1.0
QMOYAN_positive (station, horizon) float64 48B 0.689 ... 0.7263
QMOYAN_changed_positive (station, horizon) float64 48B 0.689 ... 0.7263
QMOYAN_negative (station, horizon) float64 48B 0.311 ... 0.2737
QMOYAN_changed_negative (station, horizon) float64 48B 0.311 ... 0.2737
QMOYAN_valid (station, horizon) float64 48B 0.09524 ... 0.09524
QMOYAN_agree (station, horizon) float64 48B 0.689 ... 0.7263
Attributes: (12/37)
title: Simulation Atlas Hydroclimatique 2020
summary: Derive relative des indicateurs hydrologi...
institution: DEH (Direction de l'Expertise Hydrique)
institute_id: DEH
contact: atlas.hydroclimatique@environnement.gouv....
experiment: Hydrological simulation
... ...
cat:bias_adjust_project: INFO-Crue-CM5
cat:mip_era: CMIP5
cat:experiment: rcp45
cat:bias_adjust_institution: OURANOS
cat:id: INFO-Crue-CM5_CMIP5_rcp45
ensemble_size: 378Cela entraîne une grande quantité de données avec de nombreuses variables uniques. Pour simplifier les résultats, nous allons calculer trois nouvelles statistiques :
Le changement médian.
L’écart interquartile du changement.
L’accord entre les modèles en utilisant les catégories du GIEC.
La dernière statistique est légèrement plus complexe. Pour plus de détails sur les catégories d’accord sur le signe du changement, consultez le résumé technique dans « 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.
Pour le calculer, vous pouvez utiliser les résultats produits par robustness_fractions, mais cela nécessite un appel à la fonction xclim.ensembles.robustness_categories. Les seuils et opérations requièrent deux entrées : la première est liée au test de significativité, et la seconde fait référence au pourcentage de simulations montrant un delta positif. Par exemple, « Accord envers une augmentation » est satisfait si plus de 66% des simulations montrent un changement significatif, et 80% des simulations montrent un changement positif.
[17]:
out = xr.Dataset()
out["QMOYAN_median"] = ens_stats["QMOYAN"].sel(percentiles=50)
out["QMOYAN_iqr"] = ens_stats["QMOYAN"].sel(percentiles=75) - ens_stats["QMOYAN"].sel(
percentiles=25
)
[18]:
from xclim.ensembles import robustness_categories
categories = [
"Agreement towards increase",
"Agreement towards decrease",
"Conflicting signals",
"No change or robust signal",
]
thresholds = [[0.66, 0.8], [0.66, 0.2], [0.66, 0.8], [0.66, np.nan]]
ops = [[">=", ">="], [">=", "<="], [">=", "<"], ["<", None]]
out["QMOYAN_robustness_categories"] = robustness_categories(
changed_or_fractions=ens_stats["QMOYAN_changed"],
agree=ens_stats["QMOYAN_positive"],
categories=categories,
thresholds=thresholds,
ops=ops,
)
Enfin, en utilisant une méthode de perturbation, les valeurs futures pour QMOYAN peuvent être obtenues en multipliant l’indicateur de référence avec les percentiles des delta de l’ensemble.
[19]:
out["QMOYAN_projected"] = ref.QMOYAN * (1 + ens_stats.QMOYAN / 100)
[20]:
out
[20]:
<xarray.Dataset> Size: 972B
Dimensions: (station: 2, horizon: 3, percentiles: 5)
Coordinates:
lat (station) float32 8B 49.53 49.53
lon (station) float32 8B -77.05 -77.02
drainage_area (station) float32 8B 60.76 54.61
* station (station) <U9 72B 'ABIT00057' 'ABIT00058'
* horizon (horizon) <U9 108B '2021-2051' ... '2071-2100'
* percentiles (percentiles) int64 40B 10 25 50 75 90
realization <U86 344B ...
Data variables:
QMOYAN_median (station, horizon) float64 48B 3.078 ... 4.254
QMOYAN_iqr (station, horizon) float64 48B 6.981 ... 9.842
QMOYAN_robustness_categories (station, horizon) int64 48B 3 3 3 3 3 3
QMOYAN_projected (station, percentiles, horizon) float64 240B ...9.2.3. Statistiques d’ensemble avec des données de référence probabilistes¶
Cette méthode est similaire à la section précédente, mais elle s’applique à des cas comme l”Atlas hydroclimatique du Québec méridional ou les résultats du Notebook Interpolation optimale, où les indicateurs hydrologiques pour la période historique sont représentés par une fonction de densité de probabilité (PDF) plutôt que par une seule valeur discrète. Dans de tels cas, les percentiles de l’ensemble ne peuvent pas simplement être multipliés par la valeur de référence.
Dans cet exemple, au lieu d’une seule valeur, QMOYAN est représenté par 21 percentiles qui capturent l’incertitude entourant cette statistique. Comme pour les simulations futures, nous avons également 6 plateformes hydrologiques à prendre en compte.
AVERTISSEMENT
Dans ces cas, les percentiles dans ref représentent l’incertitude (par exemple, liée à la modélisation hydrologique ou à l’incertitude des données d’entrée), et non la variabilité interannuelle. À ce stade, la moyenne temporelle devrait déjà avoir été calculée.
[21]:
ref = xr.open_mfdataset(reference_files, combine="nested", concat_dim="platform")
ref
[21]:
<xarray.Dataset> Size: 3kB
Dimensions: (platform: 6, station: 2, percentile: 21)
Coordinates:
lat (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
lon (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
* percentile (percentile) float32 84B 1.0 5.0 10.0 15.0 ... 90.0 95.0 99.0
drainage_area (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
* station (station) <U9 72B 'ABIT00057' 'ABIT00058'
realization (platform) <U86 2kB 'A20_HYDREP_QCMERI_XXX_DEBITJ_HIS_XXX_...
Dimensions without coordinates: platform
Data variables:
QMOYAN (platform, station, percentile) float32 1kB dask.array<chunksize=(1, 2, 21), meta=np.ndarray>
Attributes: (12/29)
title: Analyse des indicateurs hydrologiques aux...
summary: Hydrological indicator analysis
institution: DEH (Direction de l'Expertise Hydrique)
institute_id: DEH
contact: edouard.mailhot@environnement.gouv.qc.ca
date_created: 2020-10-30
... ...
cat:frequency: fx
cat:variable: QMOYAN
cat:type: reconstruction-hydro
cat:processing_level: indicators
cat:xrfreq: fx
cat:id: PortraitHydro_HYDROTEL_MG24HS_GovQc_GCQ_QC[22]:
# This can also be represented as a cumulative distribution function (CDF)
import matplotlib.pyplot as plt
for platform in ref.platform:
plt.plot(
ref.QMOYAN.isel(station=0).sel(platform=platform),
ref.QMOYAN.percentile / 100,
"grey",
)
plt.xlabel("Mean annual flow (m³/s)")
plt.ylabel("Probability")
plt.title("CDF for QMOYAN @ ABIT00057 \nEach line is an hydrological platform")
En raison de leur nature probabiliste, les valeurs de référence historiques ne peuvent pas être facilement combinées avec les deltas futurs. Pour résoudre ce problème, les fonctions xhydro.cc.weighted_random_sampling et xhydro.cc.sampled_indicators ont été conçues. Ensemble, ces fonctions vont :
Échantillonner “n” valeurs de la distribution historique, en fonction de la dimension “percentile”.
Échantillonner “n” valeurs de la distribution des deltas, en utilisant les poids fournis.
Créer la distribution future en appliquant les deltas échantillonnés à la distribution historique échantillonnée élément par élément.
Calculer les percentiles de la distribution future.
Tout d’abord, nous allons échantillonner dans le jeu de données de référence pour combiner les résultats des 6 plateformes hydrologiques ensemble.
[23]:
Help on function weighted_random_sampling in module xhydro.cc:
weighted_random_sampling(ds: xarray.core.dataset.Dataset, *, weights: xarray.core.dataarray.DataArray | None = None, include_dims: list[str] | None = None, n: int = 5000, seed: int | None = None) -> xarray.core.dataset.Dataset
Sample from a dataset using random sampling.
Parameters
----------
ds : xr.Dataset
Dataset to sample from.
weights : xr.DataArray, optional
Weights to use when sampling the dataset, for dimensions other than 'percentile' or 'quantile'.
See Notes for more information on special cases for 'percentile', 'quantile', 'time', and 'horizon' dimensions.
include_dims : list of str, optional
List of dimensions to include when sampling the dataset, in addition to the 'percentile' or 'quantile' dimensions and those with weights.
These dimensions will be sampled uniformly.
n : int
Number of samples to generate.
seed : int, optional
Seed to use for the random number generator.
Returns
-------
xr.Dataset
Dataset containing the 'n' samples, stacked along a 'sample' dimension.
Notes
-----
If the dataset contains a "percentile" [0, 100] or "quantile" [0, 1] dimension, the percentiles will be sampled accordingly as to
account for the uneven spacing between percentiles and maintain the distribution's shape.
Weights along the 'time' or 'horizon' dimensions are supported, but behave differently than other dimensions. They will not
be stacked alongside other dimensions in the new 'sample' dimension. Rather, a separate sampling will be done for each time/horizon,
[24]:
deltas = xclim.ensembles.create_ensemble(ds_dict_deltas)
hist_dist = xh.cc.weighted_random_sampling(
ds=ref,
include_dims=["platform"],
n=10000,
seed=0,
)
hist_dist
[24]:
<xarray.Dataset> Size: 4MB
Dimensions: (station: 2, sample: 10000)
Coordinates:
* station (station) <U9 72B 'ABIT00057' 'ABIT00058'
lat (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
lon (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
drainage_area (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
realization (sample) <U86 3MB dask.array<chunksize=(10000,), meta=np.ndarray>
* sample (sample) int64 80kB 0 1 2 3 4 5 ... 9995 9996 9997 9998 9999
Data variables:
QMOYAN (station, sample) float32 80kB dask.array<chunksize=(2, 10000), meta=np.ndarray>
Attributes: (12/29)
title: Analyse des indicateurs hydrologiques aux...
summary: Hydrological indicator analysis
institution: DEH (Direction de l'Expertise Hydrique)
institute_id: DEH
contact: edouard.mailhot@environnement.gouv.qc.ca
date_created: 2020-10-30
... ...
cat:frequency: fx
cat:variable: QMOYAN
cat:type: reconstruction-hydro
cat:processing_level: indicators
cat:xrfreq: fx
cat:id: PortraitHydro_HYDROTEL_MG24HS_GovQc_GCQ_QC[25]:
# Let's show how the historical distribution was sampled and reconstructed
def _make_cdf(ds, bins):
count, bins_count = np.histogram(ds.QMOYAN.isel(station=0), bins=bins)
pdf = count / sum(count)
return bins_count, np.cumsum(pdf)
# Barplot
plt.subplot(2, 1, 1)
uniquen = np.unique(hist_dist.QMOYAN.isel(station=0), return_counts=True)
plt.bar(uniquen[0], uniquen[1], width=0.01, color="k")
plt.ylabel("Number of instances")
plt.title("Sampling within the historical distribution")
# CDF
plt.subplot(2, 1, 2)
for i, platform in enumerate(ref.platform):
plt.plot(
ref.QMOYAN.isel(station=0).sel(platform=platform),
ref.percentile / 100,
"grey",
label="CDFs from the percentiles" if i == 0 else None,
)
bc, c = _make_cdf(hist_dist, bins=50)
plt.plot(bc[1:], c, "r", label=f"Sampled historical CDF (n={10000})", linewidth=3)
plt.ylabel("Probability")
plt.xlabel("QMOYAN (m³/s)")
plt.legend()
plt.tight_layout()
Nous pouvons faire de même pour les deltas. Puisque weights contient déjà toutes les dimensions que nous voulons échantillonner, nous n’avons pas besoin de include_dims ici.
[26]:
delta_dist = xh.cc.weighted_random_sampling(
ds=deltas,
weights=weights,
n=10000,
seed=0,
)
delta_dist
[26]:
<xarray.Dataset> Size: 320kB
Dimensions: (station: 2, horizon: 3, sample: 10000)
Coordinates:
* station (station) <U9 72B 'ABIT00057' 'ABIT00058'
* horizon (horizon) <U9 108B '2021-2051' '2041-2071' '2071-2100'
lat (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
lon (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
drainage_area (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
* sample (sample) int64 80kB 0 1 2 3 4 5 ... 9995 9996 9997 9998 9999
Data variables:
QMOYAN (station, horizon, sample) float32 240kB dask.array<chunksize=(2, 3, 10000), meta=np.ndarray>
Attributes: (12/44)
title: Simulation Atlas Hydroclimatique 2020
summary: Derive relative des indicateurs hydrologi...
institution: DEH (Direction de l'Expertise Hydrique)
institute_id: DEH
contact: atlas.hydroclimatique@environnement.gouv....
date_created: 2020-11-23
... ...
cat:bias_adjust_project: INFO-Crue-CM5
cat:mip_era: CMIP5
cat:activity: CMIP5
cat:experiment: rcp45
cat:member: r4i1p1
cat:bias_adjust_institution: OURANOS[27]:
# Then, let's show how the deltas were sampled, for the last horizon
plt.subplot(2, 1, 1)
uniquen = np.unique(delta_dist.QMOYAN.isel(station=0, horizon=-1), return_counts=True)
plt.bar(uniquen[0], uniquen[1], width=0.25, color="k")
plt.ylabel("Number of instances")
plt.title("Sampling within the historical distribution")
plt.subplot(2, 1, 2)
bc, c = _make_cdf(delta_dist, bins=100)
plt.plot(bc[1:], c, "k", label=f"Sampled deltas CDF (n={10000})", linewidth=3)
plt.ylabel("Probability")
plt.xlabel("Deltas (%)")
plt.legend()
plt.tight_layout()
Une fois que les deux distributions ont été acquises, xhydro.cc.sampled_indicators peut être utilisé pour les combiner élément par élément et reconstruire une distribution future. La distribution résultante possédera les dimensions uniques des deux jeux de données. Ici, cela signifie que nous obtenons une distribution reconstruite pour chaque horizon futur.
[28]:
help(xh.cc.sampled_indicators)
Help on function sampled_indicators in module xhydro.cc:
sampled_indicators(ds_dist: xarray.core.dataset.Dataset, deltas_dist: xarray.core.dataset.Dataset, *, delta_kind: str | dict | None = None, percentiles: xarray.core.dataarray.DataArray | xarray.core.dataset.Dataset | list[int] | numpy.ndarray | None = None) -> xarray.core.dataset.Dataset | tuple[xarray.core.dataset.Dataset, xarray.core.dataset.Dataset]
Reconstruct indicators using a perturbation approach and random sampling.
Parameters
----------
ds_dist : xr.Dataset
Dataset containing the sampled reference distribution, stacked along a 'sample' dimension.
Typically generated using 'xhydro.cc.weighted_random_sampling'.
deltas_dist : xr.Dataset
Dataset containing the sampled deltas, stacked along a 'sample' dimension.
Typically generated using 'xhydro.cc.weighted_random_sampling'.
delta_kind : str or dict, optional
Type of delta provided. Recognized values are: ['absolute', 'abs.', '+'], ['percentage', 'pct.', '%'].
If a dict is provided, it should map the variable names to their respective delta type.
If None, the variables should have a 'delta_kind' attribute.
percentiles : xr.DataArray or xr.Dataset or list or np.ndarray, optional
Percentiles to compute in the future distribution.
This will change the output of the function to a tuple containing the future distribution and the future percentiles.
If given a Dataset, it should contain a 'percentile' or 'quantile' dimension.
If given a list or np.ndarray, it should contain percentiles [0, 100].
Returns
-------
fut_dist : xr.Dataset
The future distribution, stacked along the 'sample' dimension.
fut_pct : xr.Dataset
Dataset containing the future percentiles.
Notes
-----
The future percentiles are computed as follows:
1. Sample 'n' values from the reference distribution.
2. Sample 'n' values from the delta distribution.
3. Create the future distribution by applying the sampled deltas to the sampled historical distribution, element-wise.
4. Compute the percentiles of the future distribution (optional).
[29]:
fut_dist, fut_pct = xh.cc.sampled_indicators(
ds_dist=hist_dist,
deltas_dist=delta_dist,
delta_kind="percentage",
percentiles=ref.percentile,
)
fut_dist
WARNING: utils.clean_up(): Unable to generate a new id for the dataset. Got single positional indexer is out-of-bounds.
[29]:
<xarray.Dataset> Size: 4MB
Dimensions: (station: 2, sample: 10000, horizon: 3)
Coordinates:
* station (station) <U9 72B 'ABIT00057' 'ABIT00058'
lat (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
lon (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
drainage_area (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
realization (sample) <U86 3MB dask.array<chunksize=(10000,), meta=np.ndarray>
* sample (sample) int64 80kB 0 1 2 3 4 5 ... 9995 9996 9997 9998 9999
* horizon (horizon) <U9 108B '2021-2051' '2041-2071' '2071-2100'
Data variables:
QMOYAN (station, sample, horizon) float32 240kB dask.array<chunksize=(2, 10000, 3), meta=np.ndarray>Puisque nous avons utilisé l’argument percentiles, il a également calculé une série de percentiles.
[30]:
fut_pct
[30]:
<xarray.Dataset> Size: 1kB
Dimensions: (station: 2, horizon: 3, percentile: 21)
Coordinates:
* station (station) <U9 72B 'ABIT00057' 'ABIT00058'
lat (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
lon (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
drainage_area (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
* horizon (horizon) <U9 108B '2021-2051' '2041-2071' '2071-2100'
* percentile (percentile) float32 84B 1.0 5.0 10.0 15.0 ... 90.0 95.0 99.0
Data variables:
QMOYAN (percentile, station, horizon) float64 1kB dask.array<chunksize=(21, 2, 3), meta=np.ndarray>[31]:
# The distributions themselves can be used to create boxplots and compare the historical distribution to the future ones.
plt.boxplot(
[
hist_dist.QMOYAN.isel(station=0),
fut_dist.QMOYAN.isel(station=0, horizon=0),
fut_dist.QMOYAN.isel(station=0, horizon=1),
fut_dist.QMOYAN.isel(station=0, horizon=2),
],
labels=["Historical", "2011-2040", "2041-2070", "2071-2100"],
)
plt.ylabel("Mean summer flow (m³/s)")
plt.tight_layout()
Les mêmes statistiques qu’auparavant peuvent également être calculées en utilisant les 10 000 échantillons dans delta_dist.
[32]:
# The same statistics as before can also be computed by using delta_dist
delta_dist = delta_dist.rename({"sample": "realization"}) # xclim compatibility
ens_stats = xh.cc.ensemble_stats(delta_dist, statistics)
out_prob = xr.Dataset()
out_prob["QMOYAN_median"] = ens_stats["QMOYAN"].sel(percentiles=50)
out_prob["QMOYAN_iqr"] = ens_stats["QMOYAN"].sel(percentiles=75) - ens_stats[
"QMOYAN"
].sel(percentiles=25)
out_prob["QMOYAN_robustness_categories"] = robustness_categories(
changed_or_fractions=ens_stats["QMOYAN_changed"],
agree=ens_stats["QMOYAN_positive"],
categories=categories,
thresholds=thresholds,
ops=ops,
)
out_prob
[32]:
<xarray.Dataset> Size: 308B
Dimensions: (station: 2, horizon: 3)
Coordinates:
* station (station) <U9 72B 'ABIT00057' 'ABIT00058'
* horizon (horizon) <U9 108B '2021-2051' ... '2071-2100'
lat (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
lon (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
drainage_area (station) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
percentiles int64 8B 50
Data variables:
QMOYAN_median (station, horizon) float32 24B dask.array<chunksize=(2, 3), meta=np.ndarray>
QMOYAN_iqr (station, horizon) float32 24B dask.array<chunksize=(2, 3), meta=np.ndarray>
QMOYAN_robustness_categories (station, horizon) int64 48B dask.array<chunksize=(2, 3), meta=np.ndarray>