7. Regional frequency analyses¶

[1]:

import matplotlib.pyplot as plt
import xdatasets
from lmoments3.distr import KappaGen
from sklearn.cluster import HDBSCAN, OPTICS, AgglomerativeClustering

import xhydro as xh
import xhydro.frequency_analysis as xhfa
import xhydro.gis as xhgis

ERROR 1: PROJ: proj_create_from_database: Open of /home/docs/checkouts/readthedocs.org/user_builds/xhydro/conda/stable/share/proj failed

Redefining 'percent' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining '%' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining 'year' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining 'yr' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining 'C' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining 'd' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining 'h' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining 'degrees_north' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining 'degrees_east' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining 'degrees' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining '[speed]' (<class 'pint.delegates.txt_defparser.plain.DerivedDimensionDefinition'>)

/home/docs/checkouts/readthedocs.org/user_builds/xhydro/conda/stable/lib/python3.14/site-packages/xhydro/__init__.py:21: UserWarning: The `exactextract` library is not present in the environment and will not be used.

This notebook will demonstrate how to use the xHydro library to perform regional frequency analyses on a dataset of streamflow data. Since the initial steps for regional frequency analysis are the same as those for local frequency analysis, users are encouraged to refer to the Local frequency analysis notebook for an overview.

For this example, we will use the same dataset of hydrometric gauges covering parts of southern Quebec, ensuring continuity with the previous analysis while expanding to a regional scale. However, as regional analyses rely on having access to multiple sources of data, we will extract streamflow for 15 stations.

[2]:

ds = (
    xdatasets.Query(
        **{
            "datasets": {
                "deh": {
                    "id": ["02*"],
                    "regulated": ["Natural"],
                    "variables": ["streamflow"],
                }
            },
            "time": {"start": "1970-01-01", "minimum_duration": (30 * 365, "d")},
        }
    )
    .data.squeeze()
    .load()
)

# This dataset lacks some attributes, so let's add them.
ds = ds.rename({"streamflow": "q"})
ds["id"].attrs["cf_role"] = "timeseries_id"
ds["q"].attrs = {
    "long_name": "Streamflow",
    "units": "m3 s-1",
    "standard_name": "water_volume_transport_in_river_channel",
    "cell_methods": "time: mean",
}

# Clean some of the coordinates that are not needed for this example
ds = ds.drop_vars([c for c in ds.coords if c not in ["id", "time", "name"]])

timeargs = {
    "annual": {},
}

ds_4fa = xh.indicators.get_yearly_op(
    ds, op="max", timeargs=timeargs, missing="pct", missing_options={"tolerance": 0.15}
)
ds_4fa

7.1. Explanatory variables¶

In regional frequency analyses, explanatory variables are used to help explain the spatial variation in hydrological extremes across different locations. These variables can include factors such as catchment area, elevation, precipitation, land use, and soil type, among others. By incorporating explanatory variables, we can account for the influence of geographic and environmental characteristics on extreme events, allowing for more accurate regional predictions.

7.1.2. b) Principal component analysis¶

After acquiring the explanatory variables, the next step is to process them using Principal Component Analysis (PCA) to reduce the dimensionality of the dataset. PCA helps simplify the dataset by transforming the original variables into a smaller set of uncorrelated components while retaining most of the variation in the data. This is accomplished using the function xhydro.frequency_analysis.regional.fit_pca.

[4]:

help(xhfa.regional.fit_pca)

Help on function fit_pca in module xhydro.frequency_analysis.regional:

fit_pca(ds: xr.Dataset, **kwargs) -> tuple
    Perform Principal Component Analysis (PCA) on the input dataset.

    This function scales the input data, applies PCA transformation, and returns
    the transformed data along with the PCA object.

    Parameters
    ----------
    ds : xr.Dataset
        Input dataset to perform PCA on.
    \*\*kwargs : dict
        Additional keyword arguments to pass to the PCA constructor.

    Returns
    -------
    tuple: A tuple containing:
        - data_pca (xr.DataArray): PCA-transformed data with 'Station' and 'components' as coordinates.
        - obj_pca (sklearn.decomposition.PCA): Fitted PCA object.

    Notes
    -----
    - The input data is scaled before PCA is applied.
    - The number of components in the output depends on the n_components parameter passed to PCA.

[5]:

data_pca, pca = xhfa.regional.fit_pca(dswp, n_components=3)
data_pca

[5]:

<xarray.DataArray (Station: 15, components: 3)> Size: 360B
array([[-2.48506812e+03,  8.68222390e-01,  6.18492284e-01],
       [-2.28706847e+03, -3.82682765e-01, -7.46540544e-01],
       [-1.48206850e+03,  3.89329737e-01, -2.41928874e-01],
       [-3.82066870e+02, -4.09517146e-01,  2.10528949e-01],
       [-1.85066639e+02,  2.97231987e-01,  3.20226530e-03],
       [ 1.12933711e+02, -2.73945787e+00, -2.80370822e-01],
       [ 4.13933839e+02,  1.43184979e+00,  4.50497174e-02],
       [ 5.11933876e+02,  9.70448285e-01,  1.23118865e+00],
       [ 5.32933768e+02, -5.49290384e-01, -5.67606131e-01],
       [ 5.38934065e+02, -3.95607723e-01, -8.15645783e-01],
       [ 5.42933642e+02, -1.23796579e+00,  1.14791838e+00],
       [ 8.11933845e+02, -6.70129061e-01,  8.80067651e-01],
       [ 1.11393391e+03,  2.05320024e-01, -4.18140383e-01],
       [ 1.11793429e+03,  4.39850876e+00, -5.36706342e-01],
       [ 1.12393366e+03, -2.17626024e+00, -5.29509020e-01]])
Coordinates:
  * Station     (Station) object 120B '020404' '020602' ... '024007' '024013'
  * components  (components) int64 24B 0 1 2
Attributes:
    long_name:         Fitted Scaled Data
    description:       Fitted scaled data with StandardScaler and PCA from sk...
    fitted_variables:  ['area', 'perimeter', 'gravelius', 'centroid_lon', 'ce...

The results show that the correlation between the components is close to 0, indicating that the first three components are sufficiently independent. This suggests that these components can be effectively used for the remainder of our analysis, as they capture most of the variation in the explanatory variables without significant overlap or multicollinearity.

[6]:

data_pca.to_dataframe(name="value").reset_index().pivot(
    index="Station", columns="components"
).corr()

[6]:

		value
	components	0	1	2
	components
value	0	1.000000e+00	9.942103e-14	1.350617e-14
	1	9.942103e-14	1.000000e+00	-4.653426e-16
	2	1.350617e-14	-4.653426e-16	1.000000e+00

7.1.3. c) Clustering¶

The results from the PCA can be used to group the stations into clusters based on their similarities in the principal components. Clustering helps identify regions with similar characteristics, enabling more targeted and accurate regional frequency analyses. This step is accomplished using xhydro.frequency_analysis.regional.get_clusters, which supports clustering methods from sklearn.cluster. In this example, we will use AgglomerativeClustering to form 3 clusters based on the PCA results.

[7]:

help(xhfa.regional.get_clusters)

Help on function get_clusters in module xhydro.frequency_analysis.regional:

get_clusters(model: Callable, param: dict, sample: xr.Dataset | xr.DataArray) -> list
    Get indices of groups from a fit using the specified model and parameters.

    Parameters
    ----------
    model : callable
        Model class or instance with a fit method.
    param : dict
        Parameters for the model.
    sample : xr.Dataset or xr.DataArray
        Data sample to fit the model.

    Returns
    -------
    list :
        List of indices for each non-excluded group.

[8]:

groups = xhfa.regional.get_clusters(
    AgglomerativeClustering, {"n_clusters": 3}, data_pca
)
groups

[8]:

[array(['023303', '023401', '023422', '023428', '023432', '023701',
        '024003', '024007', '024013'], dtype=object),
 array(['020404', '020602', '021407'], dtype=object),
 array(['022507', '022704', '023002'], dtype=object)]

[9]:

ax = plt.subplot(1, 1, 1)
gdf[gdf["Station"].isin(groups[0])].plot(ax=ax, color="red")
gdf[gdf["Station"].isin(groups[1])].plot(ax=ax, color="green")
gdf[gdf["Station"].isin(groups[2])].plot(ax=ax, color="blue")

[9]:

<Axes: >

../_images/notebooks_regional_frequency_analysis_14_1.png

7. Regional frequency analyses¶

7.1. Explanatory variables¶

7.1.1. a) Extraction of watershed characteristics using `xhydro.gis`¶

7.1.2. b) Principal component analysis¶

7.1.3. c) Clustering¶

7.2. Regional Frequency Analysis¶

7.3. Uncertainties¶

7.3.1. a) Bootstrapping the observations¶

7.3.2. b) Using multiple regions¶

7.3.3. c) Combining the bootstrap with multiple regions¶

7.3.4. d) Comparison¶

7. Regional frequency analyses¶

7.1. Explanatory variables¶

7.1.1. a) Extraction of watershed characteristics using xhydro.gis¶

7.1.2. b) Principal component analysis¶

7.1.3. c) Clustering¶

7.2. Regional Frequency Analysis¶

7.3. Uncertainties¶

7.3.1. a) Bootstrapping the observations¶

7.3.2. b) Using multiple regions¶

7.3.3. c) Combining the bootstrap with multiple regions¶

7.3.4. d) Comparison¶

7.1.1. a) Extraction of watershed characteristics using `xhydro.gis`¶