Sobol Sensitivity Analysis for SWAT+ Model Parameters¶

This notebook demonstrates how to perform a Sobol sensitivity analysis on SWAT+ model parameters using the pySWATPlus package and the SALib library. The analysis focuses on two parameters, epco and esco, in the hydrology.hyd file.

---

1. Import Required Libraries¶

First, we import the necessary libraries, including pySWATPlus for SWAT+ model interaction, SALib for sensitivity analysis, and concurrent.futures for parallel execution.

In [ ]:

from pySWATPlus.TxtinoutReader import TxtinoutReader
import random
from SALib.sample import saltelli
from SALib.analyze import sobol
import numpy as np
import random
import concurrent.futures

from pySWATPlus.TxtinoutReader import TxtinoutReader import random from SALib.sample import saltelli from SALib.analyze import sobol import numpy as np import random import concurrent.futures

---

2. Initialize the TxtinoutReader¶

We initialize the TxtinoutReader to interact with the SWAT+ project files.

In [5]:

txtinout_reader = TxtinoutReader('/mnt/c/Users/joans/OneDrive/Escriptori/icra/muga_windows')

txtinout_reader = TxtinoutReader('/mnt/c/Users/joans/OneDrive/Escriptori/icra/muga_windows')

---

3. Configure the SWAT+ Simulation¶

We configure the SWAT+ simulation by setting the simulation period, warmup period, and enabling the output of the channel_sd variable.

In [6]:

txtinout_reader.set_beginning_and_end_year(2010, 2012)  # Set simulation period
txtinout_reader.set_warmup(1)  # Set warmup period
txtinout_reader.enable_object_in_print_prt('channel_sd', True, False, False, False)  # Enable output for 'channel_sd'

txtinout_reader.set_beginning_and_end_year(2010, 2012) # Set simulation period txtinout_reader.set_warmup(1) # Set warmup period txtinout_reader.enable_object_in_print_prt('channel_sd', True, False, False, False) # Enable output for 'channel_sd'

---

4. Define the Model Evaluation Function¶

We define a function to run the SWAT+ model with specific values for epco and esco and evaluate the model's output. This function will be used in the sensitivity analysis.

In [7]:

def run_and_evaluate_swat(epco: float = 0.5, esco: float = 0.5):
    """
    Run the SWAT+ model with specific values for `epco` and `esco` and evaluate the output.

    Parameters:
        epco (float): Plant evaporation compensation factor.
        esco (float): Soil evaporation compensation factor.

    Returns:
        float: A mock error metric (to be replaced with actual evaluation logic).
    """
    print(f'Running SWAT with epco = {epco} and esco = {esco} \\n')

    # Copy the SWAT+ project to a temporary directory for parallel execution
    tmp_path = txtinout_reader.copy_swat()
    reader = TxtinoutReader(tmp_path)

    # Register and modify the 'hydrology.hyd' file
    hydrology_hyd = reader.register_file('hydrology.hyd', has_units=False)
    hydrology_hyd_df = hydrology_hyd.df

    # Overwrite the values of `epco` and `esco`
    hydrology_hyd_df['epco'] = epco
    hydrology_hyd_df['esco'] = esco

    # Save the changes
    hydrology_hyd.overwrite_file()

    # Run the SWAT+ model
    txt_in_out_result = reader.run_swat(show_output=False)

    # Read the results
    result_reader = TxtinoutReader(txt_in_out_result)
    channel_sdmorph = result_reader.register_file('channel_sdmorph_day.txt', has_units=True)
    channel_sdmorph_df = channel_sdmorph.df

    # Here, you should read your observations and calculate the objective function
    # For now, we return a mock value
    return random.random()

# Wrapper function for parallel execution
def evaluate(params):
    """
    Wrapper function for parallel execution of `run_and_evaluate_swat`.

    Parameters:
        params (tuple): A tuple containing the values for `epco` and `esco`.

    Returns:
        float: The result of `run_and_evaluate_swat`.
    """
    return run_and_evaluate_swat(*params)

def run_and_evaluate_swat(epco: float = 0.5, esco: float = 0.5): """ Run the SWAT+ model with specific values for `epco` and `esco` and evaluate the output. Parameters: epco (float): Plant evaporation compensation factor. esco (float): Soil evaporation compensation factor. Returns: float: A mock error metric (to be replaced with actual evaluation logic). """ print(f'Running SWAT with epco = {epco} and esco = {esco} \\n') # Copy the SWAT+ project to a temporary directory for parallel execution tmp_path = txtinout_reader.copy_swat() reader = TxtinoutReader(tmp_path) # Register and modify the 'hydrology.hyd' file hydrology_hyd = reader.register_file('hydrology.hyd', has_units=False) hydrology_hyd_df = hydrology_hyd.df # Overwrite the values of `epco` and `esco` hydrology_hyd_df['epco'] = epco hydrology_hyd_df['esco'] = esco # Save the changes hydrology_hyd.overwrite_file() # Run the SWAT+ model txt_in_out_result = reader.run_swat(show_output=False) # Read the results result_reader = TxtinoutReader(txt_in_out_result) channel_sdmorph = result_reader.register_file('channel_sdmorph_day.txt', has_units=True) channel_sdmorph_df = channel_sdmorph.df # Here, you should read your observations and calculate the objective function # For now, we return a mock value return random.random() # Wrapper function for parallel execution def evaluate(params): """ Wrapper function for parallel execution of `run_and_evaluate_swat`. Parameters: params (tuple): A tuple containing the values for `epco` and `esco`. Returns: float: The result of `run_and_evaluate_swat`. """ return run_and_evaluate_swat(*params)

---

5. Define the Sensitivity Analysis Problem¶

We define the problem for the Sobol sensitivity analysis, specifying the parameters (epco and esco) and their bounds.

In [8]:

problem = {
    'num_vars': 2,  # Number of parameters
    'names': ['epco', 'esco'],  # Parameter names
    'bounds': [[0, 1]] * 2  # Parameter bounds
}

problem = { 'num_vars': 2, # Number of parameters 'names': ['epco', 'esco'], # Parameter names 'bounds': [[0, 1]] * 2 # Parameter bounds }

---

6. Generate Parameter Samples¶

We generate parameter samples using the Saltelli sampling method from the SALib library.

In [9]:

param_values = saltelli.sample(problem, 2)  # Generate parameter samples

param_values = saltelli.sample(problem, 2) # Generate parameter samples

/tmp/ipykernel_387433/2689363164.py:1: DeprecationWarning: `salib.sample.saltelli` will be removed in SALib 1.5.1 Please use `salib.sample.sobol`
  param_values = saltelli.sample(problem, 2)  # Generate parameter samples

---

7. Perform Parallel Model Evaluations¶

We use parallel processing to evaluate the SWAT+ model for each set of parameter values.

In [10]:

# Parallel execution of model evaluations
with concurrent.futures.ProcessPoolExecutor() as executor:
    y = np.array(list(executor.map(evaluate, param_values)))

# Parallel execution of model evaluations with concurrent.futures.ProcessPoolExecutor() as executor: y = np.array(list(executor.map(evaluate, param_values)))

Running SWAT with epco = 0.46875 and esco = 0.46875 \nRunning SWAT with epco = 0.46875 and esco = 0.46875 \nRunning SWAT with epco = 0.09375 and esco = 0.65625 \nRunning SWAT with epco = 0.09375 and esco = 0.46875 \nRunning SWAT with epco = 0.59375 and esco = 0.15625 \nRunning SWAT with epco = 0.46875 and esco = 0.65625 \nRunning SWAT with epco = 0.59375 and esco = 0.96875 \nRunning SWAT with epco = 0.96875 and esco = 0.96875 \nRunning SWAT with epco = 0.96875 and esco = 0.96875 \nRunning SWAT with epco = 0.96875 and esco = 0.15625 \nRunning SWAT with epco = 0.09375 and esco = 0.65625 \n










Running SWAT with epco = 0.59375 and esco = 0.15625 \n

---

8. Analyze the Results¶

We perform the Sobol sensitivity analysis to calculate the sensitivity indices.

In [11]:

sobol_indices = sobol.analyze(problem, y)  # Perform Sobol analysis

sobol_indices = sobol.analyze(problem, y) # Perform Sobol analysis

/home/zephol/miniconda3/envs/pyswatplus/lib/python3.12/site-packages/SALib/util/__init__.py:274: FutureWarning: unique with argument that is not not a Series, Index, ExtensionArray, or np.ndarray is deprecated and will raise in a future version.
  names = list(pd.unique(groups))

---

9. Interpret the Results¶

The sobol_indices object contains the first-order, second-order, and total-order sensitivity indices for epco and esco. These indices can be used to understand the relative importance of each parameter in influencing the model output.

In [12]:

print("First-order Sobol indices:", sobol_indices['S1'])
print("Total-order Sobol indices:", sobol_indices['ST'])

print("First-order Sobol indices:", sobol_indices['S1']) print("Total-order Sobol indices:", sobol_indices['ST'])

First-order Sobol indices: [-0.84706442  1.00546731]
Total-order Sobol indices: [0.55593247 1.05855786]