3. Module d’interpolation optimale¶

L’interpolation optimale est un outil qui permet de combiner un champ spatialement distribué (c’est-à-dire le « champ de fond ») avec des observations ponctuelles de telle sorte que l’ensemble du champ puisse être ajusté en fonction des écarts entre les observations et le champ au point d’observation. Par exemple, il peut être utilisé pour combiner un champ de précipitation d’une” réanalyse (par exemple ERA5) avec des enregistrements d’observation, et ainsi ajuster les précipitations de la réanalyse sur l’ensemble du domaine de manière statistiquement optimale.

This page demonstrates how to use xhydro to perform optimal interpolation using field-like simulations and point observations for hydrological modelling. In this case, the background field is a set of outputs from a distributed hydrological model and the observations correspond to real hydrometric stations. The aim is to correct the background field (i.e. the distributed hydrological simulations) using optimal interpolation, as in Lachance-Cloutier et al. (2017).

Lachance-Cloutier, S., Turcotte, R. and Cyr, J.F., 2017. Combining streamflow observations and hydrologic simulations for the retrospective estimation of daily streamflow for ungauged rivers in southern Quebec (Canada). Journal of hydrology, 550, pp.294-306.

[1]:

import datetime as dt
from functools import partial
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pooch
import xarray as xr
from scipy.stats import norm

import xhydro as xh
import xhydro.optimal_interpolation.optimal_interpolation_fun as opt
from xhydro.optimal_interpolation.ECF_climate_correction import general_ecf
from xhydro.testing.helpers import deveraux

/home/docs/checkouts/readthedocs.org/user_builds/xhydro-fr/conda/latest/lib/python3.12/site-packages/esmpy/interface/loadESMF.py:94: VersionWarning: ESMF installation version 8.8.0, ESMPy version 8.8.0b0
  warnings.warn("ESMF installation version {}, ESMPy version {}".format(
ERROR 1: PROJ: proj_create_from_database: Open of /home/docs/checkouts/readthedocs.org/user_builds/xhydro-fr/conda/latest/share/proj failed

3.1. Un exemple rapide¶

Imaginez un scénario dans lequel nous disposons de 3 stations d’observation des débits et d’un modèle hydrologique qui simule les débits sur ces 3 sites et sur 2 autres sites supplémentaires (pour un total de 5 sites de simulation). Nous souhaitons améliorer la qualité des simulations sur chacun des 5 sites et encore plus sur les 2 sites supplémentaires où il n’y a pas d’observations pour aider à entraîner le modèle. La configuration pourrait ressembler à ceci :

Station 1 : Observé + simulé
Station 2 : Observé + simulé
Station 3 : Observé + simulé
Station 4 : Simulation uniquement
Station 5 : Simulation uniquement

L’interpolation optimale peut aider, mais nous aurons besoin de quelques informations de base concernant chacune des stations (simulées et observées) :

Taille du bassin versant (pour mettre à l’échelle les erreurs)
Latitudes/longitudes du bassin versant, pour développer le modèle d’erreur spatiale
Données observées aux 3 emplacements jaugés
Données simulées sur les 5 sites

Définissons-les maintenant et montrons les stations sur une carte :

[2]:

# Define the station coordinates. Start with observed (Stations 1-3)
lat_obs = np.array([45.0, 45.5, 45.8])
lon_obs = np.array([-71.0, -70.4, -71.8])

# And now the simulated (Stations 1-5). Notice that the first three stations are the same as
# the observation stations.
lat_est = np.array([45.0, 45.5, 45.8, 44.2, 45.4])
lon_est = np.array([-71.0, -70.4, -71.8, -70.2, -71.9])

# We need catchment areas to scale data, in the same order as the latitudes and longitudes.
drainage_area = np.array([480.0, 270.0, 620.0, 1000.0, 700.0])

# We also need streamflow data for each site. In this example, we will use a single day, but
# it would be possible to apply to each day of a time series. Notice that while the first three
# stations correspond to the same stations for both the observed_flow and simulated_flow data,
# The data is different because one is observed and the other is simulated.
observed_flow = np.array([100.0, 110.0, 150.0])
simulated_flow = np.array([108.0, 135.0, 148.0, 154.0, 88.0])

test

Nous disposons désormais des données de base nécessaires pour démarrer le traitement en utilisant une interpolation optimale. Cependant, avant de le faire, nous devons fournir quelques hyperparamètres. Certains sont plus complexes que d’autres, décomposons donc les principales étapes.

The first is the need to compute differences (also referred to as « departures ») between observations and simulations where they both occur simultaneously. We also need to scale the data by the catchment area to ensure errors are relative and can then be interpolated. We then take the logarithm of these values to ensure extrapolation does not cause negative streamflow. We will reverse the transformation later.

[3]:

# Log-transform the scaled flow data
scaled_simulated_flow = np.log(simulated_flow / drainage_area)
scaled_observed_flow = np.log(observed_flow / drainage_area[0:3])

# Compute the departure for the three observation stations.
departures = scaled_simulated_flow[0:3] - scaled_observed_flow

Nous aurons maintenant besoin de certaines informations qui pourront (ou non) être disponibles pour nos sites d’observation et de simulation. Il s’agit notamment d’estimations de :

La variance des observations sur les sites jaugés.
La variance des débits simulés aux sites d’observation. Il s’agit d’un vecteur de taille 3 dans cet exemple, soit une valeur par site d’observation. A noter que cette variance est celle des simulations sur les sites d’observation, et non la variance des observations elles-mêmes.
La variance des débits simulés aux sites d’estimation. Il s’agit d’un vecteur de taille 5 dans cet exemple, soit une valeur par point de simulation, y compris ceux qui correspondent également à un site d’observation.

Nous ne connaissons pas ces valeurs pour cet exemple de test, mais ces valeurs peuvent être estimées dans des applications réelles en utilisant de longues séries chronologiques de débits log-transformés et mis à l’échelle ou en utilisant l’erreur de mesure de l’instrumentation sur des sites jaugés. Pour cet exemple, nous supposerons des valeurs simples de 1,0 pour chaque cas.

[4]:

# Define the ratio of the observed variance to that of the simulations. We set it to 0.15, but it could be adjusted
# according to the level of confidence in each measure.
var_obs = np.array([0.15, 0.15, 0.15])

# Define the background field (simulations) variance at the observed sites
bg_var_obs = (np.array([1.0, 1.0, 1.0]),)

# Define the background field (simulations) variance at the simulated sites. Note that the first three stations
# are the same as in the previous variable, as in our test case, the first three simulated stations are located
# at the observation sites.
bg_var_est = (np.array([1.0, 1.0, 1.0, 1.0, 1.0]),)

Si nous disposions de meilleures estimations de ces variables, nous pourrions remplacer les valeurs 1,0 par des valeurs plus appropriées. Toutefois, celles-ci peuvent également être ajustées en fonction de l’expérience passée ou par essais et erreurs.

La dernière pièce du puzzle est celle de la fonction de covariance d’erreur. En un mot, l’interpolation optimale prendra en compte la distance entre une observation (ou plusieurs observations) et le site où nous devons estimer la nouvelle valeur de débit. On comprend aisément qu’une station de simulation très proche d’une station d’observation doive être fortement corrélée avec celle-ci, alors qu’un point plus éloigné serait moins corrélé. Nous avons donc besoin d’une fonction de covariance qui estime (1) le degré de covariabilité entre un point observé et simulé en fonction de (2) la distance entre ces points. Il s’agit de la fonction ECF dont plusieurs modèles existent dans la littérature. Dans de nombreux cas, une forme de modèle sera imposée et les paramètres seront ajustés de telle sorte que le modèle représente la covariance existante entre les points.

Dans cet exemple de test, nous avons trop peu de points et pas assez de pas de temps pour établir un modèle (et un paramétrage) significatif à partir des données. Nous imposons donc un modèle. Il y en a quatre qui sont intégrés à xhydro, où par[0] et par[1] sont les paramètres du modèle à calibrer (dans des circonstances normales) et où h est la distance entre les points :

Modèle 1 : par[0] * (1 + h / par[1]) * exp(-h / par[1]) – D’après Lachance-Cloutier et al. 2017.
Modèle 2 : par[0] * exp(-0,5 * (h / par[1])**2)
Modèle 3 : par[0] * exp(-h / par[1])
Modèle 4 : par[0] * exp(-(h ** par[1]) / par[0])

Nous utiliserons le modèle n°4, mais vous pouvez le modifier ci-dessous et voir comment cela affecte les résultats. Les paramètres peuvent également être modifiés pour évaluer leurs impacts.

[5]:

# Define the ECF function model. We use partial functions to do so, with the models being predefined in the
# xhydro.optimal_interpolation.ECF_climate_correction package.
ecf_fun = partial(general_ecf, form=4)

# Now we can parameterize the ecf_fun with the model parameters that we want.
# In this test example we will use values of 1.0 and 0.5 for par[0] and par[1], respectively.
par = [1.0, 0.5]
ecf = partial(ecf_fun, par=par)

Nous pouvons maintenant exécuter l’algorithme d’interpolation optimal et récupérer la valeur estimée et la variance de l’incertitude pour chaque site simulé.

[6]:

print(f"lat_est: {lat_est}")
print(f"lon_est: {lon_est}")
print(f"lat_obs: {lat_obs}")
print(f"lon_obs: {lon_obs}")
print(f"bg_departures: {departures}")
print(f"bg_est: {scaled_simulated_flow}")
print(f"bg_var_obs: {bg_var_obs}")
print(f"bg_var_est: {bg_var_est}")
print(f"var_obs: {var_obs}")
print(f"ecf: {ecf}")

lat_est: [45.  45.5 45.8 44.2 45.4]
lon_est: [-71.  -70.4 -71.8 -70.2 -71.9]
lat_obs: [45.  45.5 45.8]
lon_obs: [-71.  -70.4 -71.8]
bg_departures: [ 0.07696104  0.20479441 -0.01342302]
bg_est: [-1.49165488 -0.69314718 -1.4325072  -1.87080268 -2.07374352]
bg_var_obs: (array([1., 1., 1.]),)
bg_var_est: (array([1., 1., 1., 1., 1.]),)
var_obs: [0.15 0.15 0.15]
ecf: functools.partial(<function general_ecf at 0x71ed1d678a40>, form=4, par=[1.0, 0.5])

[7]:

# Display all the data that will be used for the optimal interpolation:


# Launch the optimal interpolation with all the pre-established values.
v_est, var_est, _ = opt.optimal_interpolation(
    lat_est=lat_est,
    lon_est=lon_est,
    lat_obs=lat_obs,
    lon_obs=lon_obs,
    bg_departures=departures,
    bg_est=scaled_simulated_flow,
    bg_var_obs=bg_var_obs,
    bg_var_est=bg_var_est,
    var_obs=var_obs,
    ecf=ecf,
    precalcs={},  # To speedup calculations, not required here.
)

Enfin, nous pouvons calculer la valeur réelle de la meilleure estimation et la variance de la distribution de l’incertitude à partir de ces résultats :

[8]:

# Transform back into absolute values and rescale by the drainage area
estimated_flow = np.exp(v_est) * drainage_area

print(f"Estimated values are: {estimated_flow}")
print(f"Simulated values were: {simulated_flow}")
print(f"Observed values are: {observed_flow}")

Estimated values are: [113.44733348 161.14678052 142.02962921 154.00037475  87.99639936]
Simulated values were: [108. 135. 148. 154.  88.]
Observed values are: [100. 110. 150.]

Dans une dernière étape, voici comment estimer la distribution des valeurs possibles aux sites d’estimation. Le v_est est l’emplacement de la distribution et le v_est est la variance. Cela signifie que nous pouvons modéliser la distribution et générer des valeurs de quantiles.

[9]:

# Get the log-normal error model, whose variance depends on the interpolation variance and the variance of the departures.
var_bg = np.var(departures)  # Variance of the departures of the background field
var_est = (
    var_est * var_bg
)  # Complete error model that includes the interpolation variance and the departure variance.

# Using the uncertainty estimation, get the 25th percentile of the estimated flows, and un-transform
percentile_values = norm.ppf(np.array(25.0) / 100.0, loc=v_est, scale=np.sqrt(var_est))
flows_25th_percentile = np.exp(percentile_values) * drainage_area

# Do the same but with the 75th percentile.
percentile_values = norm.ppf(np.array(75.0) / 100.0, loc=v_est, scale=np.sqrt(var_est))
# Get the values in real units and scale according to drainage area
flows_75th_percentile = np.exp(percentile_values) * drainage_area

print(f"Estimated values for the 25th percentile are: {flows_25th_percentile}")
print(f"Estimated values for the 50th percentile are: {estimated_flow}")
print(f"Estimated values for the 75th percentile are: {flows_75th_percentile}")

Estimated values for the 25th percentile are: [[111.26527089 158.04725979 139.29781056 144.97663776  82.84020497]]
Estimated values for the 50th percentile are: [113.44733348 161.14678052 142.02962921 154.00037475  87.99639936]
Estimated values for the 75th percentile are: [[115.67218928 164.30708704 144.81502252 163.58577347  93.473529  ]]

Notice that there are a few problems with the model presented here:

The optimal interpolation worsened the estimated results at the gauged sites compared to the raw simulation.
The 25th and 75th percentile values for the estimated flows at the gauged sites are « nan ».
The estimated flows at the ungauged sites did not change (or changed very little).

These problems arise due to some methodological choices:

Forcing of a covariance function model and parameterization that is inadequate.
Very few observation stations, making it extremely difficult to assess spatial patterns.
Simulated and observed flows that were randomly generated and have no (or little) correlation, very small covariance.

Cela signifie que le problème est mal défini et que l’interpolation optimale ne doit pas être appliquée dans ces cas. Avec plus de données, les résultats deviennent bien meilleurs, comme le montrera la section suivante.

3.2. Application sur données réelles du modèle hydrologique HYDROTEL¶

La section précédente a montré comment implémenter l’algorithme d’interpolation optimal de manière autonome. Toutefois, cela n’est pas pratique lorsque de nombreuses stations doivent être traitées simultanément. Des outils ont donc été intégrés à xhydro pour faciliter tout le traitement et, en tant que tels, ont des exigences spécifiques en matière de données. Ici, nous explorons le contenu d’un fichier d’entrée complet, et nous ajouterons des détails un peu plus tard. Commençons par importer quelques données de test depuis le dépôt xhydro-testdata :

[10]:

# Get data
test_data_path = deveraux().fetch(
    "optimal_interpolation/OI_data_corrected.zip",
    pooch.Unzip(),
)
directory_to_extract_to = Path(test_data_path[0]).parent

# Read-in all the files and set to paths that we can access later.
flow_obs_info_file = directory_to_extract_to / "A20_HYDOBS_TEST_corrected.nc"
flow_sim_info_file = directory_to_extract_to / "A20_HYDREP_TEST_corrected.nc"
corresponding_station_file = directory_to_extract_to / "station_correspondence.nc"
selected_station_file = (
    directory_to_extract_to / "stations_retenues_validation_croisee.csv"
)

We now have 4 files:

flow_obs_info_file: The dataset file (.nc) that contains the point observations and station metadata.
flow_sim_info_file: The dataset file (.nc) that contains the background field simulations, including simulated station metadata.
corresponding_station_file: The dataset file (.nc) that links the station identifiers between observations and simulated stations. This is necessary because observed stations have « real world » identifiers and distributed simulations are often coded or numbered sequentially. However, we need to be able to find which of the background field stations (simulation points) correspond to each real-world station.
selected_station_file: The list of stations from the observation set that we wish to use (thus discarding the others from the flow_obs_info_file set).

Nous pouvons maintenant les traiter pour extraire certaines valeurs qu’il faudra envoyer au contrôleur principal de l’interpolation optimale :

[11]:

# We first open the .nc files that contain our required data (qobs for the observed flows,
# qsim for the simulated flows and the station correspondence file)
qobs = xr.open_dataset(flow_obs_info_file)
qsim = xr.open_dataset(flow_sim_info_file)
station_correspondence = xr.open_dataset(corresponding_station_file)

# Also read the .csv file that contains the list of observation stations to include.
df_validation = pd.read_csv(selected_station_file, sep=None, dtype=str)
observation_stations = list(df_validation["No_station"])

Explorons le contenu de ces fichiers :

[12]:

# First show the contents of the observed streamflow file:
display(qobs)

3. Module d’interpolation optimale¶

3.1. Un exemple rapide¶

3.2. Application sur données réelles du modèle hydrologique HYDROTEL¶

3.2.1. IMPORTANT:¶

3.2.2. IMPORTANT:¶