Note
Go to the end to download the full example code.
Sensitivity analysis#
Before running a Data Assimilation it is often necessary to evaluate the sensitivity of the model parameters with respect to a given scenario. Here we use the Weil et al dataset and generate 24 trajectories varying PERMX and POROS (hydraulic conductivity and porosity of the soil).
Estimated time to run the notebook = 5min
import multiprocessing
import os
import shutil
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from SALib.analyze import morris as ma
from SALib.plotting import morris as mp
from SALib.sample import morris as ms
from pyCATHY import cathy_tools
from pyCATHY.cathy_tools import subprocess_run_multi
prj_name = "test0"
path2prj = "weil_exemple_sensitivityAnalysis" # add your local path here
simu = cathy_tools.CATHY(dirName=path2prj, prj_name=prj_name)
simu.run_preprocessor(verbose=False)
simu.run_processor(
IPRT1=2,
DTMIN=1e-2,
DTMAX=1e2,
DELTAT=5,
TRAFLAG=0,
verbose=False,
)
SPP_map = simu.set_SOIL_defaults(SPP_map_default=True)
dsw, _ = simu.read_outputs("sw")
obs_data = dsw.to_numpy() # (n_times, n_nodes)
obs = obs_data[-1, :] # last time step, (n_nodes,)
n_nodes = obs.size
🏁 Initiate CATHY object
🍳 gfortran compilation
👟 Run preprocessor
🔄 Update parm file
🔄 Update hap.in file
🔄 Update dem_parameters file
🔄 Update dem_parameters file
🛠 Recompile src files [4s]
🍳 gfortran compilation [11s]
✅ Compilation successful!
👟 Run processor
morris_problem = {
"num_vars": 2,
"names": ["PERMX", "POROS"],
"bounds": [
[1e-5, 1e-3], # Ks
[0.4, 0.6], # porosity
],
"groups": None,
}
Morris gives N * (num_vars + 1) samples: 8 * (2 + 1) = 24
number_of_trajectories = 8
sample = ms.sample(morris_problem, number_of_trajectories, num_levels=4)
df_sample = pd.DataFrame(sample, columns=morris_problem["names"])
df_sample.index.name = "sample"
# SPP_map is (zone, layer)-indexed, dtype object: cast to float and
# collapse to one reference value (homogeneous soil)
for p in morris_problem["names"]:
ref = SPP_map[p].astype(float).mean()
df_sample["dev_" + p] = 1e2 * (df_sample[p] - ref) / ref
fig = plt.figure()
mp.sample_histograms(fig, sample, morris_problem)
[(array([9., 2., 9., 4.]), array([1.000e-05, 2.575e-04, 5.050e-04, 7.525e-04, 1.000e-03]), <BarContainer object of 4 artists>), (array([5., 5., 7., 7.]), array([0.4 , 0.45, 0.5 , 0.55, 0.6 ]), <BarContainer object of 4 artists>)]
pathexe_list = []
for ii in range(len(sample)):
path_exe = os.path.join(
simu.workdir, prj_name + "_sensitivity", "sample" + str(ii + 1)
)
pathexe_list.append(path_exe)
if os.path.exists(path_exe):
continue
shutil.copytree(os.path.join(simu.workdir, prj_name), path_exe)
simu.update_veg_map()
for ii in range(len(sample)):
# fresh defaults each iteration (avoids carry-over between samples)
SoilPhysProp = simu.set_SOIL_defaults(SPP_map_default=True)
SoilPhysProp[["PERMX", "PERMY", "PERMZ"]] = sample[ii, 0]
SoilPhysProp["POROS"] = sample[ii, 1]
simu.update_soil(SPP_map=SoilPhysProp, path=pathexe_list[ii] + "/input/")
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
🔄 Update soil
homogeneous soil
with multiprocessing.Pool(processes=multiprocessing.cpu_count()) as pool:
result = pool.map(subprocess_run_multi, pathexe_list)
simu_ensemble = np.zeros((n_nodes, len(sample)))
for ii in range(len(sample)):
dsw_ii, _ = simu.read_outputs("sw", path=pathexe_list[ii] + "/output/")
simu_ensemble[:, ii] = dsw_ii.to_numpy()[-1, :]
noise = 0.025 * obs # assume 2.5% noise in the data
rmse = np.array([
err_weighted_rmse(simu_ensemble[:, ii], obs, noise)
for ii in range(len(sample))
])
assert rmse.shape == (len(sample),)
print("rmse:", rmse)
print("max spread across samples:", np.ptp(simu_ensemble, axis=1).max())
print("return codes:", result)
# 1) did the soil files actually differ between samples?
for p in pathexe_list[:3]:
print(p, os.listdir(os.path.join(p, "input")))
print(open(os.path.join(p, "input", "soil")).read()[:300]) # adjust filename if needed
# 2) are the outputs different files / freshly written?
for p in pathexe_list[:3]:
f = os.path.join(p, "output", "sw")
print(f, os.path.exists(f), os.path.getmtime(f) if os.path.exists(f) else None)