Note
Go to the end to download the full example code.
DA with Random Initial Conditions on Soil Layers#
This notebook demonstrates how to create an ensemble of models with varying initial conditions for different soil layers. By introducing randomness into the initial conditions, we can better simulate the natural variability in soil properties and improve the robustness of data assimilation (DA) processes.
Estimated time to run the notebook = 2min
import os
import numpy as np
import pyvista as pv
from pyCATHY.DA.perturbate import perturbate_parm
from pyCATHY.DA import perturbate
from pyCATHY.DA.cathy_DA import DA
from pyCATHY.plotters import cathy_plots as CTp
🏁 Initiate CATHY object
😟 src files not found
working directory
is:/home/runner/work/pycathy_wrapper/pycathy_wrapper/examples/DA/./DA_with_non_u
niform_ic
📥 Fetch cathy src files
📥 Fetch cathy prepro src files
📥 Fetch cathy input files
🔄 Update hap.in file
🔄 Update dem_parameters file
🔄 Update dem_parameters file
simuWithDA.update_dem_parameters()
simuWithDA.update_veg_map()
simuWithDA.update_soil()
🔄 Update dem_parameters file
────────────────────────── ⚠ warning messages above ⚠ ──────────────────────────
The parm dictionnary is empty
Falling back to defaults to update CATHYH
This can have consequences !!
────────────────────────────────────────────────────────────────────────────────
🔄 Update parm file
────────────────────────── ⚠ warning messages above ⚠ ──────────────────────────
['Adjusting TMAX with respect to time of interests requested\n']
────────────────────────────────────────────────────────────────────────────────
🔄 Update soil
Inconsistent number of zones/layers with respect to the number of soil lines: 15/6!Revise dem_parm file!
homogeneous soil
simuWithDA.NENS = 3
list_pert = perturbate.perturbate(
simuWithDA,
scenario,
simuWithDA.NENS,
)
var_per_dict_stacked = {}
for dp in list_pert:
var_per_dict_stacked = perturbate_parm(
var_per_dict_stacked,
parm=dp,
type_parm = dp['type_parm'], # can also be VAN GENUCHTEN PARAMETERS
mean = dp['mean'],
sd = dp['sd'],
sampling_type = dp['sampling_type'],
ensemble_size = dp['ensemble_size'], # size of the ensemble
per_type= dp['per_type'],
nlayers = nlay,
savefig= os.path.join(simuWithDA.workdir,
simuWithDA.project_name,
simuWithDA.project_name + dp['savefig'])
)
/opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/scipy/stats/_qmc.py:958: UserWarning: The balance properties of Sobol' points require n to be a power of 2.
sample = self._random(n, workers=workers)
This in normally directly called when using run_DA_sequential()
simuWithDA._create_subfolders_ensemble()
simuWithDA.update_ENS_files(var_per_dict_stacked,
var_per_dict_stacked.keys(),
cycle_nb=0
)
😟 processor exe not found
🔄 Update ensemble
🍳 gfortran compilation
👟 Run preprocessor
🔄 Update parm file
🛠 Recompile src files [7s]
🍳 gfortran compilation [13s]
b''
👟 Run processor
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
Single value detected for ic ==> assumming it homogeneous
pl = pv.Plotter(shape=(1,2))
for i, ensi in enumerate([1,3]):
DApath = f'DA_Ensemble/cathy_{ensi}/vtk/'
path = os.path.join(simuWithDA.workdir,
simuWithDA.project_name,
DApath,
simuWithDA.project_name + '.vtk'
)
pl.subplot(0,i)
CTp.show_vtk(path,
'ic',
ax=pl,
clim = [-25,-5],
show_scalar_bar=True,
)
_ = pl.add_legend('')
pl.add_title(f'Ensemble nb:{ensi}')
pl.show()
plot ic
plot ic
Total running time of the script: (0 minutes 15.799 seconds)